DataStructs.h

Go to the documentation of this file.

//
//
// TCAlg data structs
//
// Bruce Baller
//
#ifndef TRAJCLUSTERALGDATASTRUCT_H
#define TRAJCLUSTERALGDATASTRUCT_H

// C/C++ standard libraries
#include <array>
#include <bitset>
#include <vector>

// LArSoft libraries
#include "larcorealg/Geometry/fwd.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"

namespace TMVA {
  class Reader;
}
namespace calo {
  class CalorimetryAlg;
}

namespace tca {

  namespace detail {
    struct SortEntry {
      unsigned int index;
      float val;
    };
  }

  using Point3_t = std::array<double, 3>;
  using Vector3_t = std::array<double, 3>;
  using Point2_t = std::array<float, 2>;
  using Vector2_t = std::array<double, 2>;

  // some functions to handle the CTP_t type
  typedef unsigned int CTP_t;
  constexpr unsigned int Tpad = 10;    // alignment for CTP sub-items - TPC
  constexpr unsigned int Cpad = 10000; // alignment for CTP sub-items - Cryostat

  inline CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
  {
    return cryo * Cpad + tpc * Tpad + plane;
  }
  inline CTP_t EncodeCTP(const geo::PlaneID& planeID)
  {
    return EncodeCTP(planeID.Cryostat, planeID.TPC, planeID.Plane);
  }
  inline CTP_t EncodeCTP(const geo::WireID& wireID)
  {
    return EncodeCTP(wireID.Cryostat, wireID.TPC, wireID.Plane);
  }
  geo::PlaneID DecodeCTP(CTP_t CTP);


  struct VtxStore {
    Point2_t Pos{{0, 0}};
    Point2_t PosErr{{2, 1}};
    unsigned short NTraj{0};
    unsigned short Pass{0}; // Pass in which this vertex was created
    float ChiDOF{0};
    // Topo: 0 = end0-end0, 1 = end0(1)-end1(0), 2 = end1-end1, 3 = CI3DV,
    //       4 = C3DIVIG, 5 = FHV, 6 = FHV2, 7 = SHCH, 8 = CTBC, 9 = Junk, 10 = HamBragg, 11 = neutral decay (pizero)
    //       12 = BraggSplit, 13 = Not Used, 14 = Reconcile2VTs
    short Topo{0};
    CTP_t CTP{0};
    int ID{0};  
    int UID{0}; 
    int Vx3ID{0};
    float Score{0};
    float TjChgFrac{0}; 
    std::bitset<16> Stat{0}; 
  };

  typedef enum {
    kVtxTrjTried, 
    kFixed,       
    kOnDeadWire,
    kHiVx3Score, 
    kVxTruMatch, 
    kVxMerged,
    kVxIndPlnNoChg, 
    kVxEnvOK,       
    kVxBitSize      
  } VtxBit_t;

  struct Vtx3Store {
    float X{0};      // x position
    float XErr{0.5}; // x position error
    float Y{0};      // y position
    float YErr{0.5}; // y position error
    float Z{0};      // z position
    float ZErr{0.5}; // z position error
    float Score{0};
    int Wire{-1}; // wire number for an incomplete 3D vertex
    geo::TPCID TPCID;
    std::vector<int> Vx2ID; // List of 2D vertex IDs in each plane
    int ID{0};              // 0 = obsolete vertex
    int UID{0};             
    bool Primary{false};
    bool Neutrino{false};
  };

  // The (y,z) position of a valid wire intersection of two wires on two planes +
  // + the (y,z) position variation for a one wire change in pln1 and pln2
  struct TCWireIntersection {
    unsigned short pln1;
    unsigned short pln2;
    float wir1;
    float wir2;
    double y;
    double z;
    double dydw1;
    double dzdw1;
    double dydw2;
    double dzdw2;
    unsigned int tpc;
  };

  // A temporary struct for matching trajectory points; 1 struct for each TP for
  // each trajectory. These are put into mallTraj which is then sorted by increasing xlo
  struct Tj2Pt {
    unsigned int wire;
    // x range spanned by hits on the TP
    float xlo;
    float xhi;
    unsigned short plane;
    // the Trajectory ID
    unsigned short id;
    unsigned short ipt; // The trajectory point
    // the number of points in the Tj so that the minimum Tj length cut (MatchCuts[2]) can be made
    unsigned short npts;
  };

  struct TrajPoint {
    CTP_t CTP{0};            
    Point2_t HitPos{{0, 0}}; // Charge weighted position of hits in wire equivalent units
    Point2_t Pos{{0, 0}};    // Trajectory position in wire equivalent units
    Vector2_t Dir{{0, 0}};   // Direction cosines in the StepDir direction
    double HitPosErr2{0};    // Uncertainty^2 of the hit position perpendiclar to the direction
    // HitPosErr2 < 0 = HitPos not defined because no hits used
    double Ang{0};                  // Trajectory angle (-pi, +pi)
    double AngErr{0.1};             // Trajectory angle error
    float KinkSig{-1};              // kink significance = delta angle / angle error
    float Chg{0};                   // Chargetj2pt
    float AveChg{-1};               // Average charge of last ~20 TPs
    float ChgPull{0.1};             //  = (Chg - AveChg) / ChgRMS
    float Delta{0};                 // Deviation between trajectory and hits (WSE)
    float DeltaRMS{0.02};           // RMS of Deviation between trajectory and hits (WSE)
    float FitChi{0};                // Chi/DOF of the fit
    int InPFP{0};                   // ID of the PFParticle that owns this TP
    unsigned short NTPsFit{2};      // Number of trajectory points fitted to make this point
    unsigned short Step{0};         // Step number at which this TP was created
    unsigned short AngleCode{0};    // 0 = small angle, 1 = large angle, 2 = very large angle
    std::vector<unsigned int> Hits; // vector of fHits indices
    std::bitset<16> UseHit{0};      // set true if the hit is used in the fit
    std::bitset<8> Environment{0};  // TPEnvironment_t bitset that describes the environment
  };

  // struct filled by FitPar which fits TP Chg, Delta, etc
  struct ParFit {
    Point2_t Pos{{0, 0}}; // position origin of the fit
    float Par0;
    float AvePar;
    float ParErr;
    float ParSlp; // slope relative to the origin
    float ParSlpErr;
    float ChiDOF;
    unsigned short nPtsFit;
  };

  // Global information for the trajectory
  struct Trajectory {
    std::vector<TrajPoint> Pts; 
    CTP_t CTP{0};               
    std::bitset<128> AlgMod;    
    int WorkID{0};
    int ParentID{
      -1}; 
    float AveChg{0};       
    float TotChg{0};       
    float ChgRMS{0.5};     
    short MCSMom{0};       //< Crude 2D estimate to use for shower-like vs track-like discrimination
    Point2_t dEdx{{0, 0}}; 
    std::array<unsigned short, 2> VtxID{{0, 0}}; 
    std::array<unsigned short, 2> EndPt{
      {0, 0}};                 
    int ID;                    
    int UID;                   
    int SSID{0};               
    unsigned short PDGCode{0}; 
    unsigned short Pass{0};    
    short StepDir{0};          
    short StartEnd{-1};        
    std::array<std::bitset<8>, 2> EndFlag{}; // Bitset that encodes the reason for stopping
    std::bitset<8> Strategy{};               
    bool NeedsUpdate{false};                 
    bool IsGood{true}; 
    bool MaskedLastTP{false};
  };

  struct TjForecast {
    unsigned short nextForecastUpdate{
      0}; 
    float showerLikeFraction{
      0};              
    float outlook{-1}; 
    float chgSlope{0};
    float chgSlopeErr{0};
    float chgFitChiDOF{0};
    float chgRMS{0};
    short MCSMom{0};
    bool leavesBeforeEnd{false}; 
    bool foundShower{false};
    bool endBraggPeak{false};
  };

  // Temporary 3D trajectory points composed of triplet or doublet wire intersections
  struct TrajPoint3 {
    Point3_t Pos{{0.0, 0.0, 0.0}};
    Vector3_t Dir{{0.0, 0.0, 0.0}};
    std::vector<Tj2Pt> Tj2Pts; // list of trajectory points
  };

  struct SectionFit {
    Point3_t Pos{{-10.0, 0.0, 0.0}};   
    Vector3_t Dir{{0.0, 0.0, 0.0}};    
    Vector3_t DirErr{{0.0, 0.0, 0.0}}; 
    float ChiDOF{-1};
    unsigned short NPts{0};
    bool NeedsUpdate{true}; 
  };

  // a 3D trajectory point composed of a 3D point & direction and a single TP
  struct TP3D {
    Point3_t Pos{{-10.0, -10.0, -10.0}}; 
    Vector3_t Dir{{0.0, 0.0, 0.0}};
    double TPX{-10};   
    double TPXErr2{1}; 
    float Wire{-1};
    float along{1E6}; 
    int TjID{0};      
    CTP_t CTP;
    unsigned short TPIndex{USHRT_MAX}; 
    unsigned short SFIndex{USHRT_MAX}; 
    std::bitset<8> Flags;              //< see TP3DFlags_t
  };

  typedef enum {
    kTP3DGood,  // Is good for fitting and calorimetry
    kTP3DBad,   // Should be removed from the trajectory
    kTP3DHiDEdx // Has high dE/dx
  } TP3DFlags_t;

  // Struct for 3D trajectory matching
  struct MatchStruct {
    // IDs of Trajectories that match in all planes
    std::vector<int> TjIDs;
    // Count of the number of X-matched hits and de-weight by angle
    float Count{0}; // Set to 0 if matching failed
  };

  constexpr unsigned int pAlgModSize = 6; // alignment for CTP sub-items - TPC

  struct PFPStruct {
    std::vector<int> TjIDs;  // used to reference Tjs within a slice
    std::vector<int> TjUIDs; // used to reference Tjs in any slice
    std::vector<TP3D> TP3Ds; // vector of 3D trajectory points
    std::vector<SectionFit> SectionFits;
    // Start is 0, End is 1
    std::array<std::vector<float>, 2> dEdx;
    std::array<std::vector<float>, 2> dEdxErr;
    std::array<int, 2> Vx3ID{{0, 0}};
    int BestPlane{-1};
    int PDGCode{-1};
    std::vector<int> DtrUIDs;
    size_t ParentUID{0}; // Parent PFP UID (or 0 if no parent exists)
    geo::TPCID TPCID;
    float CosmicScore{0};
    int ID{0};
    int UID{0};                              // unique global ID
    unsigned short MVI;                      // matVec index for detailed debugging
    std::bitset<8> Flags;                    //< see PFPFlags_t
    std::array<std::bitset<8>, 2> EndFlag{}; // Uses the same enum as Trajectory EndFlag
    std::bitset<pAlgModSize> AlgMod; //< Allocate the first set of bits in AlgBit_t for 3D algs
  };

  typedef enum { kCanSection, kNeedsUpdate, kSmallAngle } PFPFlags_t;

  struct ShowerPoint {
    Point2_t Pos; // Hit Position in the normal coordinate system
    Point2_t
      RotPos;     // Position rotated into the shower coordinate system (0 = along, 1 = transverse)
    float Chg{0}; // Charge of this point
    unsigned int HitIndex; // the hit index
    unsigned short
      TID; // The ID of the tj the point (hit) is in. TODO eliminate this redundant variable
  };

  // A temporary structure that defines a 2D shower-like cluster of trajectories
  struct ShowerStruct {
    CTP_t CTP;
    int ShowerTjID{0};              // ID of the shower Trajectory composed of many InShower Tjs
    std::vector<int> TjIDs;         // list of InShower Tjs
    std::vector<int> NearTjIDs;     // list of Tjs that are not InShower but satisfy the maxSep cut
    std::vector<ShowerPoint> ShPts; // Trajectory points inside the shower
    float Angle{0};                 // Angle of the shower axis
    float AngleErr{3};              // Error
    float AspectRatio{1}; // The ratio of charge weighted transverse/longitudinal positions
    float DirectionFOM{1};
    std::vector<Point2_t> Envelope; // Vertices of a polygon that encompasses the shower
    float EnvelopeArea{0};
    float ChgDensity{0}; // Charge density inside the Envelope
    float Energy{0};
    float ParentFOM{10};
    int ID{0};
    int UID{0};      
    int ParentID{0}; // The ID of a parent Tj - the one at the start of the shower
    int TruParentID{0};
    int SS3ID{0};           // ID of a ShowerStruct3D to which this 2D shower is matched
    bool NeedsUpdate{true}; // Needs to be updated (e.g. after adding a tj, defining a parent, etc)
  };

  // Shower variables filled in MakeShowers. These are in cm and radians
  struct ShowerStruct3D {
    Vector3_t Dir; //
    Vector3_t
      DirErr; // DirErr is hijacked to store the shower rms at the start, center and end sections
    Point3_t Start;    //
    Point3_t StartErr; // PosErr is hijacked to temporarily store the charge in the three sections
    Point3_t ChgPos;   // position of the center of charge
    Point3_t End;      // end position
    double Len{1};
    double OpenAngle{0.12};
    std::vector<double> Energy;
    std::vector<double> EnergyErr;
    std::vector<double> MIPEnergy;
    std::vector<double> MIPEnergyErr;
    std::vector<double> dEdx;
    std::vector<double> dEdxErr;
    geo::TPCID TPCID;
    std::vector<int> CotIDs; // list of indices of 2D showers in tjs.cots
    std::vector<unsigned int> Hits;
    int BestPlane;
    int ID;
    int UID{0};      
    int ParentID{0}; // The ID of a track-like pfp at the start of the shower, e.g. an electron
    float MatchFOM;
    unsigned short PFPIndex{USHRT_MAX}; // The index of the pfp for this shower
    int Vx3ID{0};
    bool NeedsUpdate{true}; // This is set true whenever the shower needs to be updated
    bool Cheat{false};
  };

  struct DontClusterStruct {
    std::array<int, 2>
      TjIDs;   // pairs of Tjs that shouldn't be clustered in shower reconstruction because...
    int Vx2ID; // they share a 2D vertex that may be matched to...
    int Vx3ID; // a high-score 3D vertex
  };

  struct ShowerTreeVars {
    // run, subrun, and event are also saved to this tree

    std::vector<float> BeginWir; // begin wire
    std::vector<float> BeginTim; // begin tick
    std::vector<float> BeginAng; // begin angle
    std::vector<float> BeginChg; // beginning average charge
    std::vector<short> BeginVtx; // ID of begin vertex
    std::vector<float> EndWir;   // end wire
    std::vector<float> EndTim;   // end tick
    std::vector<float> EndAng;   // end angle
    std::vector<float> EndChg;   // ending average charge
    std::vector<short> EndVtx;   //ID of end vertex

    std::vector<short> MCSMom;

    std::vector<short> PlaneNum;

    std::vector<int> TjID;
    std::vector<int> IsShowerTj;        // indicates tj is an shower trajectory
    std::vector<int> ShowerID;          // shower ID associated w/ trajectory. -1 = no shower
    std::vector<int> IsShowerParent;    // this tj was chosen as a parent tj
    std::vector<int> StageNum;          // stage of reconstruction
    std::vector<std::string> StageName; // stage name

    // envelope information
    std::vector<float> Envelope;
    std::vector<int> EnvPlane;
    std::vector<int> EnvStage;
    std::vector<int> EnvShowerID;

    int nStages{0};
    unsigned short nPlanes{0};
  };

  struct CRTreeVars {
    std::vector<int> cr_origin;
    std::vector<float> cr_pfpxmin;
    std::vector<float> cr_pfpxmax;
    std::vector<float> cr_pfpyzmindis;
  };

  // Algorithm modification bits
  typedef enum {
    kFillGaps3D, // 3D algorithms for PFPs (bitset size limited to 8 bits)
    kKink3D,
    kTEP3D,
    kJunk3D,
    kRTPs3D,
    kMat3D, // 2D algorithms for Tjs here and below
    kMaskHits,
    kMaskBadTPs,
    kMichel,
    kDeltaRay,
    kCTStepChk,
    kRvPrp,
    kCHMUH,
    kSplit,
    kComp3DVx,
    kComp3DVxIG,
    kHamBragg,
    kHamVx,
    kHamVx2,
    kJunkVx,
    kJunkTj,
    kKilled,
    kMerge,
    kLastEndMerge,
    kTEP,
    kTHCEP,
    kEndKink,
    kCHMEH,
    kFillGaps,
    kUseGhostHits,
    kMrgGhost,
    kChkInTraj,
    kStopBadFits,
    kFixBegin,
    kFTBChg,
    kBeginChg,
    kFixEnd,
    kBraggSplit,
    kUUH,
    kVtxTj,
    kChkVxTj,
    kPhoton,
    kHaloTj,
    kNoFitToVx,
    kVxMerge,
    kVxNeutral,
    kNoKinkChk,
    kChkStop,
    kChkStopEP,
    kChkChgAsym,
    kFTBRvProp,
    kTjHiVx3Score,
    kVtxHitsSwap,
    kSplitHiChgHits,
    kShowerLike,
    kKillInShowerVx,
    kShowerTj,
    kShwrParent,
    kMergeOverlap,
    kMergeSubShowers,
    kMergeSubShowersTj,
    kMergeNrShowers,
    kMergeShChain,
    kCompleteShower,
    kSplitTjCVx,
    kMakePFPTjs,
    kStopShort,
    kReconcile2Vs,
    kFTBMod,
    kAlgBitSize 
  } AlgBit_t;

  typedef enum {
    kNormal,
    kStiffEl, 
    kStiffMu, 
    kSlowing  
  } Strategy_t;

  // Stop flag bits
  typedef enum {
    kSignal,
    kAtKink,
    kAtVtx,
    kBragg,
    kAtTj,
    kOutFV,
    kNoFitVx,
    kFlagBitSize 
  } EndFlag_t;

  // Environment near a trajectory point
  typedef enum {
    kEnvNotGoodWire,
    kEnvNearMuon,
    kEnvNearShower,
    kEnvOverlap,
    kEnvUnusedHits,
    kEnvNearSrcHit, 
    kEnvFlag 
  } TPEnvironment_t;

  // TrajClusterAlg configuration bits
  typedef enum {
    kStepDir,    
    kTestBeam,   
    kDebug,      
    kStudy1,     
    kStudy2,     
    kStudy3,     
    kStudy4,     
    kSaveCRTree, 
    kTagCosmics, 
    kSaveShowerTree 
  } TCModes_t;

  extern const std::vector<std::string> AlgBitNames;
  extern const std::vector<std::string> EndFlagNames;
  extern const std::vector<std::string> VtxBitNames;
  extern const std::vector<std::string> StrategyBitNames;

  // struct for configuration - used in all slices
  struct TCConfig {
    std::vector<float> vtx2DCuts; 
    std::vector<float> vtx3DCuts; 
    std::vector<float> vtxScoreWeights;
    std::vector<float> neutralVxCuts;
    std::vector<short> deltaRayTag;
    std::vector<short> muonTag;
    std::vector<float> electronTag;
    std::vector<float> chkStopCuts; 
    std::vector<float> showerTag;   
    std::vector<float> kinkCuts;    
    std::vector<float> match3DCuts; 
                                    //    std::vector<float> matchTruth;     ///< Match to MC truth
    std::vector<float> chargeCuts;
    std::vector<float> qualityCuts;           
    std::vector<float> pfpStitchCuts;         
    std::vector<unsigned short> minPtsFit;    
    std::vector<unsigned short> minPts;       
    std::vector<unsigned short> maxAngleCode; 
    std::vector<short> minMCSMom;             
    std::vector<float> angleRanges;           
    float wirePitch;
    float unitsPerTick; 
    std::vector<float> maxPos0;
    std::vector<float> maxPos1;
    float multHitSep; 
    float maxChi;
    const geo::GeometryCore* geom;
    const geo::WireReadoutGeom* wireReadoutGeom;
    calo::CalorimetryAlg* caloAlg;
    TMVA::Reader* showerParentReader;
    std::vector<float> showerParentVars;
    float hitErrFac;
    float maxWireSkipNoSignal;   
    float maxWireSkipWithSignal; 
    float projectionErrFactor;
    float VLAStepSize;
    float JTMaxHitSep2;      
    std::bitset<128> useAlg; 
    std::bitset<128> dbgAlg; 
    short recoSlice{0};      
    short recoTPC{-1};       
    bool dbgSlc{true};       
    bool dbgStp{false};      
    bool dbgMrg{false};
    bool dbg2V{false}; 
    bool dbgVxNeutral{false};
    bool dbgVxMerge{false};
    bool dbgVxJunk{false};
    bool dbg3V{false}; 
    bool dbgPFP{false};
    bool dbgDeltaRayTag{false};
    bool dbgMuonTag{false};
    bool dbg2S{false};
    bool dbg3S{false};
    bool dbgStitch{false};  
    bool dbgSummary{false}; 
    bool dbgDump{false};    
    short nPtsAve;          
    std::bitset<16> modes;  
    bool doForecast{true};
    bool useChannelStatus{true};
  };

  struct TCHit {
    unsigned int allHitsIndex; // index into fHits
    int InTraj{
      0}; // ID of the trajectory this hit is used in, 0 = none, < 0 = Tj under construction
  };

  // hit collection for all slices, TPCs and cryostats + event information
  // Note: Ideally this hit collection would be the FULL hit collection before cosmic removal
  struct TCEvent {
    geo::TPCID TPCID;
    std::vector<recob::Hit> const* allHits = nullptr;
    std::vector<recob::Hit> const* srcHits = nullptr;
    //  hit Range for the allHits collection in the current TPCID
    //   plane      wire             firstHit, lastHit
    std::vector<std::vector<std::pair<unsigned int, unsigned int>>> wireHitRange;
    // hit range for the srcHits collection
    std::vector<std::pair<unsigned int, unsigned int>> tpcSrcHitRange;
    // list of good wires in the current TPCID
    std::vector<std::vector<bool>> goodWire;
    std::vector<recob::SpacePoint> const* sptHandle =
      nullptr;                                        
    std::vector<std::array<unsigned int, 3>> sptHits; 
    unsigned int event;
    unsigned int run;
    unsigned int subRun;
    unsigned int eventsProcessed;
    std::vector<float> aveHitRMS; 
    std::vector<TCWireIntersection> wireIntersections;
    int WorkID;
    int globalT_UID;
    int globalP_UID;
    int global2V_UID;
    int global3V_UID;
    int global2S_UID;
    int global3S_UID;
    bool aveHitRMSValid{false}; 
    bool expectSlicedHits{
      false}; 
  };

  struct TCSlice {
    std::vector<unsigned int> nWires;
    std::vector<unsigned int> firstWire; 
    std::vector<unsigned int> lastWire;  
    float xLo;                           // fiducial volume of the current tpc
    float xHi;
    float yLo;
    float yHi;
    float zLo;
    float zHi;
    geo::TPCID TPCID;
    unsigned short nPlanes;
    int ID; 
    // The variables below do change in size from event to event

    // Save histograms to develop cosmic removal tools
    CRTreeVars crt;
    std::vector<TCHit> slHits;
    std::vector<Trajectory> tjs; 
    std::vector<Tj2Pt> mallTraj; 
    // vector of pairs of first (.first) and last+1 (.second) hit on each wire
    // in the range fFirstWire to fLastWire. A value of UINT_MAX indicates that there
    // are no hits on the wire.
    std::vector<std::vector<std::pair<unsigned int, unsigned int>>> wireHitRange;
    std::vector<VtxStore> vtxs;   
    std::vector<Vtx3Store> vtx3s; 
    std::vector<PFPStruct> pfps;
    std::vector<ShowerStruct> cots; // Clusters of Trajectories that define 2D showers
    std::vector<DontClusterStruct>
      dontCluster;                       // pairs of Tjs that shouldn't clustered in one shower
    std::vector<ShowerStruct3D> showers; // 3D showers
    bool isValid{false};                 // set false if this slice failed reconstruction
  };

  extern TCEvent evt;
  extern TCConfig tcc;
  extern ShowerTreeVars stv;
  extern std::vector<TjForecast> tjfs;

  // vector of hits, tjs, etc in each slice
  extern std::vector<TCSlice> slices;
  // vector of seed TPs
  extern std::vector<TrajPoint> seeds;

} // namespace tca

#endif // ifndef TRAJCLUSTERALGDATASTRUCT_H