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 <vector>
#include <bitset>

// LArSoft libraries
#include "larcore/Geometry/Geometry.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Shower.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "nusimdata/SimulationBase/MCParticle.h"

#include "TMVA/Reader.h"

namespace tca {
  
  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 = 3D split, 11 = neutral decay (pizero)
    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,      
    kVtxTruMatch,      
    kVtxMerged,
    kVtxIndPlnNoChg,  
    kVtxBitSize     
  } 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};
    short 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};
  };
  
  // 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{
    Vector2_t dir;
    unsigned int wire;
    // x range spanned by hits on the TP
    float xlo;
    float xhi;
    CTP_t ctp;
    // 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;
    short score; // 0 = Tj with nice vertex, 1 = high quality Tj, 2 = normal, -1 = already matched
  };

  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 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
    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, e.g. nearby showers or other Tjs
  };
  
  // 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};                 
    float DirFOM {0.5};         
    short MCSMom {0};         //< Crude 2D estimate to use for shower-like vs track-like discrimination
    float EffPur {0};                     
    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};               
    unsigned int mcpListIndex {UINT_MAX};
    std::array<std::bitset<8>, 2> StopFlag {};  // 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};
  };
  
  // struct used for TrajCluster 3D trajectory points
  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
    float dEdx {0};             // The charge is stored here before dE/dx is calculated
    Point2_t AlongTrans;         // Longitudinal & transverse position (cm) relative to the trajectory
  };

  // Struct for 3D trajectory matching
  struct MatchStruct {
    // IDs of Trajectories that match in all planes
    std::vector<int> TjIDs;
    std::vector<float> TjCompleteness;  // fraction of TP points that are 3D-matched
    // Count of the number of X-matched hits and de-weight by angle
    float Count {0};                    // Set to 0 if matching failed
    Point3_t Pos;               // Position center using 3D-matched points on the Tjs - 3D fit
    Vector3_t Dir;              // Direction using 3D-matched points on the Tjs - 3D fit
  };
  
  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<float> TjCompleteness;  // fraction of TP points that are 3D-matched
    std::vector<TrajPoint3> Tp3s;    // TrajCluster 3D trajectory points
    // Start is 0, End is 1
    std::array<Point3_t, 2> XYZ;        // XYZ position at both ends (cm)
    std::array<Vector3_t, 2> Dir;
    std::array<Vector3_t, 2> DirErr;
    std::array<std::vector<float>, 2> dEdx;
    std::array<std::vector<float>, 2> dEdxErr;
    std::array<int, 2> Vx3ID {{ 0, 0 }};
    int BestPlane {-1};
    // stuff for constructing the PFParticle
    int PDGCode {-1};
    std::vector<int> DtrUIDs;
    size_t ParentUID {0};       // Parent PFP UID (or 0 if no parent exists)
    geo::TPCID TPCID;
    float EffPur {0};                     
    unsigned int mcpListIndex {UINT_MAX};
    unsigned short MatchVecIndex {USHRT_MAX};
    float CosmicScore{0};
    int ID {0};
    int UID {0};              // unique global ID
    bool Primary;             // PFParticle is attached to a primary vertex
    bool NeedsUpdate {true};    // Set true if the PFParticle needs to be (re-)defined
  };

  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 {
    kMaskHits,
    kMaskBadTPs,
    kMichel,
    kDeltaRay,
    kCTKink,        
    kCTStepChk,
    kTryWithNextPass,
    kRvPrp,
    kCHMUH,
    kSplit,
    kComp3DVx,
    kComp3DVxIG,
    kHED, // High End Delta
    kHamVx,
    kHamVx2,
    kJunkVx,
    kJunkTj,
    kKilled,
    kMerge,
    kTEP,
    kCHMEH,
    kFillGap,
    kUseGhostHits,
    kMrgGhost,
    kChkInTraj,
    kStopBadFits,
    kFixBegin,
    kFTBChg,
    kBeginChg,
    kFixEnd,
    kUUH,
    kVtxTj,
    kChkVxTj,
    kMisdVxTj,
    kPhoton,
    kHaloTj,
    kNoFitToVx,
    kVxMerge,
    kVxNeutral,
    kNoKinkChk,
    kSoftKink,
    kChkStop,
    kChkStopEP,
    kChkChgAsym,
    kFTBRvProp,
    kStopAtTj,
    kMat3D,
    kM3DVxTj,
    kMat3DMerge,
    kSplit3DKink,
    kTjHiVx3Score,
    kVtxHitsSwap,
    kSplitHiChgHits,
    kShowerLike,
    kKillInShowerVx,
    kShowerTj,
    kShwrParent,
    kMergeOverlap,
    kMergeSubShowers,
    kMergeSubShowersTj,
    kMergeNrShowers,
    kMergeShChain,
    kCompleteShower,
    kSplitTjCVx,
    kNewStpCuts,
    kNewVtxCuts,
    kAlgBitSize     
  } AlgBit_t;
  
  typedef enum {
    kNormal,
    kStiffEl,       
    kStiffMu,       
    kSlowing        
  } Strategy_t;
  
  // Stop flag bits
  typedef enum {
    kSignal,
    kAtKink,
    kAtVtx,
    kBragg,
    kAtTj,
    kOutFV,
    kFlagBitSize     
  } StopFlag_t; 
  
  // Environment near a trajectory point
  typedef enum {
    kEnvDeadWire,
    kEnvNearTj,
    kEnvNearShower,
    kEnvOverlap,
    kEnvUnusedHits,
    kEnvClean,      
    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> StopFlagNames;
  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;     
    std::vector<float> chargeCuts;
    std::vector<float> qualityCuts; 
    std::vector<float> pfpStitchCuts;      
    std::vector<float> testBeamCuts;
    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 detinfo::DetectorProperties* detprop;
    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};     
    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 {false};
  };

  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 {
    std::vector<recob::Hit> const* allHits = nullptr;
    unsigned int event;
    unsigned int run;
    unsigned int subRun;
    unsigned int eventsProcessed;
    std::vector<float> aveHitRMS;      
    int WorkID;
    int globalTjID;
    int globalPFPID;
    int globalVx2ID;
    int globalVx3ID;
    int globalS2ID;
    int globalS3ID;
    bool aveHitRMSValid {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 -2 indicates that there
    // are no hits on the wire. A value of -1 indicates that the wire is dead
    std::vector<std::vector< std::pair<int, int>>> wireHitRange;
    std::vector< VtxStore > vtxs; 
    std::vector< Vtx3Store > vtx3s; 
    std::vector<MatchStruct> matchVec; 
    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