TTruthTrajectoriesModule.hxx

Go to the documentation of this file.

///
/// For questions or suggestions about this module please contact the
/// current responsible and CC in the eventAnalysis package manager.
///
/// 26-Nov-2012: Current responsible for this module is,
/// Steve Dennis (stephen.dennis [*a*t*] warwick.ac.uk)
///

#ifndef TTruthTrajectoriesModule_hxx_seen
#define TTruthTrajectoriesModule_hxx_seen

// c++
#include <set>

// Root
#include <TClonesArray.h>
#include <TLorentzVector.h>

// oaEvent
#include <TG4Trajectory.hxx>
#include <TND280Event.hxx>

// analysisMasterTools
#include "eventAnalysisEnums.hxx"

// eventAnalysis
#include <EeventAnalysis.hxx>
#include <TAnalysisTruthModuleBase.hxx>

namespace ND {
class TTruthTrajectoriesModule;
OA_EXCEPTION(EeventAnalysisInfiniteLoop, EeventAnalysis);
}

///\brief The eventAnalysis module responsible for saving information about the
/// True Particles in an event from MC.
///\details It saves a tree called "Trajectories" in the "TruthDir" directory
/// of eventAnalysis files. The tree contains integers indicating how many
/// trajectories were in the event, and a TClonesArray "Trajectories" which
/// contains a TTruthTrajectory for every particle in the event which passed the
/// module's criteria to be saved.
class ND::TTruthTrajectoriesModule : public TAnalysisTruthModuleBase {
 public:
  class TTruthTrajectory;
  class TTruthTrajectoryPoint;

  TTruthTrajectoriesModule(const char *name = "Trajectories",
                           const char *title = "True Trajectory information");

  virtual ~TTruthTrajectoriesModule();

  /// This method returns true - the module is always enabled by defult.
  virtual Bool_t IsEnabledByDefault() const { return true; }

  /// Called for the first event
  /// This method checks whether this event is a real-data event an if so throws
  /// an ND::EDataFile(). Otherwise it returns true.
  virtual Bool_t ProcessFirstEvent(ND::TND280Event &);

  /// Method for setting behaviour of module.
  /// Currently implemented
  /// options are:
  /// "saveall" - Save every trajectory longer than the minimum length
  /// regardless of where it occurred
  virtual Bool_t Configure(std::string &option);

  /// Getter for trajectory length save condition.
  Double_t GetMinimumTrajectoryLengthToSave() const { return fMinLength; }

  /// Getter for maximum number of trajectories saveable per event.
  Int_t GetMaximumNumberTrajectoriesPerEvent() const {
    return fMaxNTrajectories;
  }

  /// Getter for trajectory length save condition. [mm]
  void SetMinimumTrajectoryLengthToSave(Double_t mm) { fMinLength = mm; }

  /// Setter for maximum number of trajectories saveable per event.
  void SetMaximumNumberTrajectoriesPerEvent(UInt_t n) { fMaxNTrajectories = n; }

  ///\brief Controls whether a trajectory that would ordinarily be saved
  /// will also save all its parent trajectories.
  void SetSaveParentChain(const Bool_t &yesorno = true) {
    fSaveParentChain = yesorno;
  }

  /// Controls whether just P0D/Tracker/ECal trajectories should be saved.
  void SetSaveOnlyP0DTrackerECALTrajectories(const Bool_t &yesorno = true) {
    fSaveOnlyP0DTrackerECALTrajectories = yesorno;
  }

  ///\brief Creates the necessary tree and branches for saving the Truth
  /// Trajectories information
  virtual void InitializeBranches();

  ///\brief Called for each event, this method is the master method for
  /// retrieving and filling the Truth Trajectories information.
  virtual bool FillTree(ND::TND280Event &);

 private:
  ///\brief Fills a std::set (fSaveList) with the ID of every trajectory which
  /// should be saved for the current event.
  void FillSaveList(ND::THandle<ND::TG4TrajectoryContainer> trajectories);

  /// \brief Fills the vector of trajectory points.
  /// \details The criteria used should be controlled in the parameters file.
  void FillTrajectoryPoints(
      ND::TG4Trajectory *const traj,
      ND::TTruthTrajectoriesModule::TTruthTrajectory *trajToFill);

  ///\brief Fills vectors of entry/exit positions and momenta of the
  /// trajectory for every subdetector traversed.
  void FillTraces(ND::TG4Trajectory *const traj,
                  ND::TTruthTrajectoriesModule::TTruthTrajectory *trajToFill);

  ///\brief Determines if a TG4TrajectoryPoint is in the active region of the
  /// specified subdetector.
  ///\details Only implemented properly for TPC/FGD.
  bool GetIsActive(const ND::TG4TrajectoryPoint &point,
                   eventAnalysisEnums::ESubdetector det) const;

  /// Determines the EParticleCategory to which a trajectory belongs.
  eventAnalysisEnums::EParticleCategory GetCategory(
      const ND::TG4Trajectory *const traj);

  /// Returns true if a trajectory needs to be saved, and false oterwise.
  bool SaveTraj(ND::TG4Trajectory *const traj) const;

  ///\brief The maximum number of trajectories that can be saved from a single
  /// event.
  ///\details Initialised by the constructor.
  UInt_t fMaxNTrajectories;

  ///\brief Minimum Length of Trajectories that will be saved in mm. All primary
  /// particles will be saved regardless of this.
  Double_t fMinLength;

  ///\brief Whether to save all trajectories, or only those which intersect the
  /// P0D, Tracker or ECals.
  Bool_t fSaveOnlyP0DTrackerECALTrajectories;

  ///\brief Whether saving a trajectory should also trigger the saving of all
  /// the trajectories in its parent chain.
  Bool_t fSaveParentChain;

  /// List of the trajectory IDs which are to be saved from the current event.
  std::set<Int_t> fSaveList;

 public:
  /// Total number of trajectories saved from the event
  Int_t fNTraj;

  /// Number of Charged Lepton trajectories saved from the event
  Int_t fNTrajLepton;

  /// Number of Charged Baryon Trajectories saved from the event
  Int_t fNTrajBaryon;

  /// Number of Charged Meson Trajectories saved from the event
  Int_t fNTrajMeson;

  /// Number of Photon Trajectories saved from the event
  Int_t fNTrajPhoton;

  /// Number of Other Charged Trajectories saved from the event
  Int_t fNTrajOtherCharged;  ///< [branch]

  /// Number of Other Neutral Trajectories saved from the event
  Int_t fNTrajOtherNeutral;

  /// Number of Any Other Trajectories saved from the event
  Int_t fNTrajOther;

  ///\brief Clones array of ND::TTruthTrajectoriesModule::TTruthTrajectory
  /// sorted in ascending ID order.
  TClonesArray *fTrajectories;
};

///\brief Contains the truth information for points along the particle
/// trajectories generated during Monte Carlo simulations.
class ND::TTruthTrajectoriesModule::TTruthTrajectoryPoint : public TObject {
 public:
 
  Int_t ProcessType; ///<  The process type that caused the Monte Carlo to generate this  point.
  ///<\details The process types are detailed in the GEANT documentation, and
  ///< this field contains a packed version of the GEANT type and
  ///< sub-type. The saved value is [type + 1000*sub-type], and can be
  ///< recovered by type = ProcessType%1000 and subType = ProcessType/1000.
  ///< For quick reference, some of the interesting process type and sub-type
  ///< definitions are below.
  ///<
  ///<  * 0: The type for an undefined process.
  ///<  * 1: The type for a "transportation" process like crossing a boundary.
  ///<  * 2: The type for a electromagnetic process.
  ///<      * 1 [raw value of 1002]: The sub-type for Coulomb scattering.
  ///<      * 2 [raw value of 2002]: The sub-type for an ionization step.
  ///<      * 3 [raw value of 3002]: The sub-type for Bremsstrahlung.
  ///<      * 4 [raw value of 4002]: The sub-type for pair production by a
  ///<                charged particle.
  ///<      * 10 [raw value of 10002]: The sub-type for multiple scattering.
  ///<      * 12 [raw value of 12002]: The sub-type for the photo-electric
  ///<                effect.
  ///<      * 13 [raw value of 13002]: The sub-type for Compton scattering.
  ///<      * 14 [raw value of 14002]: The sub-type for pair production.
  ///<  * 4: The type for a hadronic process.
  ///<      * 111 [raw value of 111004]: The sub-type for hadron elastic.
  ///<      * 121 [raw value of 121004]: The sub-type for hadron inelastic.
  ///<      * 131 [raw value of 131004]: The sub-type for hadron capture.
  ///<      * 141 [raw value of 141004]: The sub-type for hadron charge
  ///<                exchange.
  ///<  * 6: The type for a decay.
  ///<  * 7: A general process (not a physics process).
  ///<      * 401 [raw value of 401007]: Step caused by a step limit.

  Float_t PositionX;   ///< The X position of the trajectory point [mm].

  Float_t PositionY;   ///< The Y position of the trajectory point [mm].

  Float_t PositionZ;   ///< The Z position of the trajectory point [mm].

  Float_t PositionT;   ///< The time of the trajectory point [ns].


  Float_t MomentumX;   ///< The X momentum of the particle leaving the trajectory point [MeV].

  Float_t MomentumY;   ///< The Y momentum of the particle leaving the trajectory point [MeV].

  Float_t MomentumZ;   ///< The Z position of the particle leaving the trajectory point [MeV].


 private:
  ClassDef(ND::TTruthTrajectoriesModule::TTruthTrajectoryPoint, 1);
};

///\brief Contains the truth information associated with a particle from
/// Monte Carlo simulations.
class ND::TTruthTrajectoriesModule::TTruthTrajectory : public TObject {
 public:
  TTruthTrajectory();
  virtual ~TTruthTrajectory(){};

  ///\brief Comparison function of trajectory IDs so that a TClonesArray can
  /// be sorted in ascending ID order.
  Int_t Compare(const TObject *obj) const;

  /// Returns true - flagging the object as being sortable.
  Bool_t IsSortable() const { return true; }

  ///\brief Trajectory's ID number. Uniquely identifies this trajectory
  /// within the event. Used by other eventAnalysis modules to reference
  /// trajectories.
  Int_t ID;

  /// Particle PDG code
  Int_t PDG;

  /// Particle name
  std::string Name;

  /// Classifier of the particle type
  Int_t Category;

  /// The Initial momentum of the particle at creation [MeV].
  TLorentzVector InitMomentum;

  /// Initial Position at which the particle was created [mm].
  TLorentzVector InitPosition;

  /// Final Position at which the particle stopped or left the simulation [mm].
  TLorentzVector FinalPosition;

  /// Parent particle's Trajectory ID.
  /// If this is a primary trajectory, the ParentID = 0.
  Int_t ParentID;

  ///\brief ID of the primary particle at the top of this particle's parent
  /// chain.
  ///\details If this is a primary trajectory, then PrimaryID = ID.
  Int_t PrimaryID;

  ///\brief Mass of the particle [MeV].
  double Mass;

  /// Charge in units of |e|/3. (e.g. an electron has charge -3)
  Int_t Charge;

  /// \brief The trajectory points where interactions occured along the
  /// trajectory.
  std::vector<ND::TTruthTrajectoriesModule::TTruthTrajectoryPoint> Points;///< The trajectory points where interactions occured along the  trajectory.

    std::vector<Int_t> TraceSubdetectors;///<\brief Vector of integer subdetector codes indicating which subdetectors the particle travelled through.
  ///<\details Uses the convention from eventAnalysisEnums::ESubdetector.
  ///< This can be an active or inactive part of that subdetector.
  ///< TraceSubdetectors.front() is the starting subdetector of the track
  ///< (i.e. equivalent to the old InitSubdetector).
  ///< TraceSubdetectors.back() is the final subdetector of the track
  ///< (i.e. equivalent to the old FinalSubdetector).
  ///< Note that all "Trace" prefixed variables store information in the
  ///< order the particle encounters the subdetector, and can deal with
  ///< encountering the same subdetector twice.


    std::vector<bool> TraceInActive;///<\brief Vector of booleans indicating whether the particle went through an active part of the corresponding subdetector.
  ///<\details TraceInActive.size() = TraceSubdetectors.size() (e.g.
  ///< TraceInActive[4] tells you whether the particle went through an active
  ///< part of the 5th subdetector in its path, TraceSubdetectors[4]).


    std::vector<TLorentzVector> TraceEntrancePosition;///<\brief Vector of TLorentzVectors that stores the entrance position of the particle in each subdetector it encounters.
  ///<\details TraceEntrancePosition.front() is the global starting position of
  ///< the track (i.e. equivalent to the old InitPosition)
  ///< TraceEntrancePosition.size() = TraceSubdetectors.size()
  ///< (e.g. TraceEntrancePosition[4] tells you the entrance position of the
  ///< 5th subdetector in its path, TraceSubdetectors[4]).


    std::vector<TLorentzVector> TraceExitPosition;///<\brief Vector of TLorentzVectors that stores the exit position of the particle in each subdetector it encounters.
  ///<\details TraceExitPosition.back() is the global stopping position of the
  ///< track (i.e. equivalent to the old FinalPosition)
  ///< TraceExitPosition.size() = TraceSubdetectors.size() (e.g.
  ///< TraceExitPosition[4] tells you the exit position of the 5th subdetector
  ///< in its path, TraceSubdetectors[4]).


    std::vector<TVector3> TraceEntranceMomentum;///<\brief Vector of TVector3s that stores the entrance momentum of the particle in each subdetector it encounters.
  ///<\details TraceEntranceMomentum.front() is the global starting momentum of
  ///< the track (i.e. equivalent to the old InitMomentum).
  ///< TraceEntranceMomentum.size() = TraceSubdetectors.size() (e.g.
  ///< TraceEntranceMomentum[4] tells you the entrance momentum of the 5th
  ///< subdetector in its path, TraceSubdetectors[4]).


    std::vector<TVector3> TraceExitMomentum;///<\brief Vector of TVector3s that stores the exit momentum of the particle in each subdetector it encounters.
  ///<\details TraceExitMomentum.back() is the global stopping momentum of the
  ///< track (i.e. equivalent to the old FinalMomentum)
  ///< TraceExitMomentum.size() = TraceSubdetectors.size() (e.g.
  ///< TraceExitMomentum[4] tells you the exit momentum of the 5th subdetector
  ///< in its path, TraceSubdetectors[4]).


 private:
  ClassDef(ND::TTruthTrajectoriesModule::TTruthTrajectory, 2);
};

#endif