DTrackCandidate_factory_CDC.h

Go to the documentation of this file.

// $Id$
//
//    File: DTrackCandidate_factory_CDC.h
// Created: Thu Sep  6 14:47:48 EDT 2007
// Creator: davidl (on Darwin Amelia.local 8.10.1 i386)
//

#ifndef _DTrackCandidate_factory_CDC_
#define _DTrackCandidate_factory_CDC_

#include <map>
#include <vector>
using namespace std;

#include "TDirectory.h"

#include <JANA/JFactory.h>
using namespace jana;

#include "DTrackCandidate.h"
#include "DHelicalFit.h"
#include "CDC/DCDCTrackHit.h"
#include <DVector3.h>
#include <HDGEOMETRY/DGeometry.h>
#include <HDGEOMETRY/DMagneticFieldMap.h>

class DTrackCandidate_factory_CDC : public JFactory<DTrackCandidate>
{
   public:
      DTrackCandidate_factory_CDC(){};
      ~DTrackCandidate_factory_CDC();
      const char* Tag(void){return "CDC";}

      enum trk_flags_t
      {
         NONE              = 0x000,
         NOISE             = 0x001,
         USED              = 0x002, //set if used in a super layer seed (whether that seed was rejected later or not)
         OUT_OF_TIME       = 0x004
      };

      enum wire_direction_t
      {
         WIRE_DIRECTION_AXIAL = 0, 
         WIRE_DIRECTION_STEREOLEFT = 1, //rotated by 6 degrees counter-clockwise
         WIRE_DIRECTION_STEREORIGHT = 2 //rotated by 6 degrees clockwise
      };

      class DCDCTrkHit
      {
         public:
            void Reset(void);
            inline double Dist2(DCDCTrkHit* trkhit) const{return (trkhit->hit->wire->origin - this->hit->wire->origin).Mag2();}

            const DCDCTrackHit* hit;
            unsigned int index;
            unsigned int flags;
            float var_z;
            DVector3 dStereoHitPos;
            float dPhiStereo;
            bool dValidStereoHitPosFlag; //false if prior to calc, or if hit doesn't intersect circle
      };

      // Collection of adjacent DCDCTrkHit's on a ring. 
      class DCDCRingSeed
      {
         public:
            //hits are stored in order from the first straw to the last straw
               //generally, this in order from lowest straw # to highest straw #
               //unless it crosses straw = 1 barrier: then first straw is lowest-# straw with phi > pi, and last straw is highest-# straw with phi < pi
            vector<DCDCTrkHit*> hits;
            bool linked;
            signed char ring;
      };

      // Collection of information about a given potential spiral turn
      typedef struct
      {
         signed char dSpiralTurnRingFlag; //-1 if potential turn looks like the track is turning back outwards, 1 if turning back inwards
         signed char dDefiniteSpiralTurnRingFlag; //same as dSpiralTurnRingFlag if DEFINITE turn, 0 otherwise
         signed char dSpiralTurnRing; //is ring# of turn
      } DSpiralParams_t;

      // Collection of adjacent DCDCRingSeed's in a super layer
      class DCDCSuperLayerSeed
      {
         public:
            void Reset(void);
            bool Are_AllHitsOnRingShared(const DCDCSuperLayerSeed* locCDCSuperLayerSeed, int locRing) const;
            inline void Get_Hits(vector<DCDCTrkHit*>& locHits) const
            {
               locHits.clear();
               for(size_t loc_i = 0; loc_i < dCDCRingSeeds.size(); ++loc_i)
                  locHits.insert(locHits.end(), dCDCRingSeeds[loc_i].hits.begin(), dCDCRingSeeds[loc_i].hits.end());
            }

            vector<DCDCRingSeed> dCDCRingSeeds; //stored in order from innermost (lowest) ring to outermost (highest) ring
            unsigned char dSuperLayer;
            unsigned int dSeedIndex;
            bool linked;
            wire_direction_t dWireOrientation;
            map<int, DSpiralParams_t> dSpiralLinkParams; //key is the dSeedIndex of the DCDCSuperLayerSeed in this super layer it is linked to (can point to itself (self-linked)!)
      };

      // Collection of adjacent DCDCSuperLayerSeed's in the CDC
         // Each unique combination of axial DCDCSuperLayerSeed's has it's own DCDCTrackCircle object
         // Every possible combination of stereo DCDCSuperLayerSeed's used to link these axial DCDCSuperLayerSeed's together is stored in the 2D-vectors
            // This is until the best combination is found: then the unused combinations are cleared and only one combination will remain
      class DCDCTrackCircle
      {
         public:
            void Reset(void);
            DCDCSuperLayerSeed* Get_LastSuperLayerSeed(void) const;
            DCDCSuperLayerSeed* Get_SuperLayerSeed(unsigned int locSuperLayer) const;
            void Strip_StereoSuperLayerSeed(unsigned int locSuperLayer);
            void Add_LastSuperLayerSeed(DCDCSuperLayerSeed* locSuperLayerSeed);
            void Truncate_Circle(unsigned int locNewLastSuperLayer);
            void Absorb_TrackCircle(const DCDCTrackCircle* locTrackCircle);
            bool Check_IfInputIsSubset(const DCDCTrackCircle* locTrackCircle);
            void Get_AllStereoSuperLayerSeeds(vector<DCDCSuperLayerSeed*>& locStereoSuperLayerSeeds);
            unsigned int Get_NumStereoSuperLayerSeeds(void);

            vector<DCDCSuperLayerSeed*> dSuperLayerSeeds_Axial;
            //for dSuperLayerSeeds_InnerStereo and dSuperLayerSeeds_OuterStereo:
               //each vector<DCDCSuperLayerSeed*> is a series of adjacent stereo seeds that could potentially belong together
               //for example: 
                  // dSuperLayerSeeds_InnerStereo[0] is likely (not necessarily) a vector with size = 2: 
                     //the first in this vector is a DCDCSuperLayerSeed from super layer 2 (SL2), the second is from SL3 and is definitely adjacent to the one from SL2.
                  // dSuperLayerSeeds_InnerStereo[1] is also likely a vector with size = 2, but a different combination of adjacent super layer seeds than the previous one. 
                  // once the best seeds have been selected, dSuperLayerSeeds_InnerStereo will have size = 1. 
               //the reason that InnerStereo and OuterStereo are kept separately is to keep track of all possible combinations of adjacent super layers
                  //that way they are evaluated together when determining which stereo super layer seeds are best
               //Note: if super layer 4 is missing (e.g. dead high voltage board), but inner (and/or outer) super layers are present, all super layers will be "inner"
                  //this is so that the stereo super layer combinatorics are evaluated correctly in the Find_ThetaZ() function
               //Note: if super layers 1 through 3 are missing (e.g. a decay product), then dSuperLayerSeeds_InnerStereo will be empty
            vector<vector<DCDCSuperLayerSeed*> > dSuperLayerSeeds_InnerStereo;
            vector<vector<DCDCSuperLayerSeed*> > dSuperLayerSeeds_OuterStereo;

            DHelicalFit* fit; //the circle fit
            float dWeightedChiSqPerDF; //of circle fit //weighted: is (chisq/ndf)/(#axial_super_layers^2): prefer fits with more axial super layers (not necessarily more hits)
            float dWeightedChiSqPerDF_Stereo; //of theta/z determination //weighted: is (chisq/ndf)/(#axial_super_layers^2): prefer fits with more stereo super layers (not necessarily more hits)
            float dAverageDriftTime; //of axial hits close to the circle fit
            vector<unsigned int> HitBitPattern; //bit pattern of hits in the track circle (first is just axial, then is all)
            float dTheta;
            float dVertexZ;
            signed char dSpiralTurnRing; //is ring# of confirmed match spiral turn, -1 if no turn

            //dTruncationSourceCircles: if this object's (not this member's) circle was truncated, this member points to the track circles that indicated truncation was necessary
               //more than 1 possible if two circles were merged after truncation:
                  //e.g. two tracks have SL1 & SL4 but different SL7, and their SL7's are rejected by other tracks who have the same SL7s
            vector<const DCDCTrackCircle*> dTruncationSourceCircles;
            bool dHasNonTruncatedSeedsFlag_InnerStereo; //true if any inner stereo seeds are present and are unique: not from a truncated source
            bool dHasNonTruncatedSeedsFlag_OuterStereo; //true if any outer stereo seeds are present and are unique: not from a truncated source
      };

      class DCDCLineFit
      {
         public:
            void BracketMinimumChisq(double &a,double &b, double &c, double &chi2a, double &chi2b, double &chi2c);
            double FindMinimumChisq(double a, double b, double c, double &lambda);
            double ChiXY(double lambda);

            unsigned int n;
            vector<float> s;
            vector<float> var_s;
            vector<float> z;
            vector<float> var_z;
            vector<float> w;
            float z0;
            float tanl;
      };

      typedef struct
      {
         float x, y, perp2;
         float z, var_z;
      } intersection_t;

      typedef vector<vector<DCDCRingSeed> >::iterator ringiter;

   private:
      jerror_t init(void);                ///< Called once at program start.
      jerror_t brun(JEventLoop *locEventLoop, int32_t runnumber); ///< Called everytime a new run number is detected.
      jerror_t evnt(JEventLoop *locEventLoop, uint64_t eventnumber); ///< Called every event.
      jerror_t erun(void);                ///< Called everytime run number changes, provided brun has been called.
      jerror_t fini(void);                ///< Called after last event of last event source has been processed.

      // Utility Functions
      void Reset_Pools(void);
      DCDCTrkHit* Get_Resource_CDCTrkHit(void);
      DCDCSuperLayerSeed* Get_Resource_CDCSuperLayerSeed(void);
      DHelicalFit* Get_Resource_HelicalFit(void);
      DCDCTrackCircle* Get_Resource_CDCTrackCircle(void);

      // Make Super Layer Seeds
      jerror_t Get_CDCHits(JEventLoop* loop);
      void Find_SuperLayerSeeds(vector<DCDCTrkHit*>& locSuperLayerHits, unsigned int locSuperLayer);
      void Link_RingSeeds(vector<DCDCRingSeed*>& parent, ringiter ring, ringiter ringend, unsigned int locSuperLayer, unsigned int locNumPreviousRingsWithoutHit);
      double MinDist2(const DCDCRingSeed& locInnerRingSeed, const DCDCRingSeed& locOuterRingSeed);
      double MinDist2(const vector<DCDCTrkHit*>& locInnerSeedHits, const vector<DCDCTrkHit*>& locOuterSeedHits);
      void Reject_SuperLayerSeeds_HighSeedDensity(unsigned int locSuperLayer);
      void Calc_SuperLayerPhiRange(DCDCSuperLayerSeed* locSuperLayerSeed, double& locSeedFirstPhi, double& locSeedLastPhi);
      bool Check_IfPhiRangesOverlap(double locFirstSeedPhi, double locLastSeedPhi, double locTargetFirstPhi, double locTargetLastPhi);

      // Search for spirals
      void Set_SpiralLinkParams(void);
      bool SearchFor_SpiralTurn_TwoSeedsSharingManyHits(DCDCSuperLayerSeed* locSuperLayerSeed1, DCDCSuperLayerSeed* locSuperLayerSeed2);
      bool SearchFor_SpiralTurn_TwoSeedsSharingFewHits(DCDCSuperLayerSeed* locSuperLayerSeed1, DCDCSuperLayerSeed* locSuperLayerSeed2);
      bool SearchFor_SpiralTurn_ManyHitsAdjacentRing(DCDCSuperLayerSeed* locSuperLayerSeed1, DCDCSuperLayerSeed* locSuperLayerSeed2, int locRingToCheck, int locMinStrawsAdjacentRing, int& locMaxSpiralNumHits);
      bool SearchFor_SpiralTurn_MissingOrBetweenRings(DCDCSuperLayerSeed* locSuperLayerSeed1, DCDCSuperLayerSeed* locSuperLayerSeed2);
      bool SearchFor_SpiralTurn_SingleSeed(DCDCSuperLayerSeed* locSuperLayerSeed);
      void Print_SuperLayerSeeds(void);

      //Link Super Layers to Create DCDCTrackCircle Objects
      bool Build_TrackCircles(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      bool Link_SuperLayers(vector<DCDCTrackCircle*>& locCDCTrackCircles, unsigned int locOuterSuperLayer);
      void Link_SuperLayers_FromAxial(vector<DCDCTrackCircle*>& locCDCTrackCircles, unsigned int locOuterSuperLayer, unsigned int locInnerSuperLayer);
      void Link_SuperLayers_FromStereo(vector<DCDCTrackCircle*>& locCDCTrackCircles, unsigned int locOuterSuperLayer, unsigned int locInnerSuperLayer);
      bool Link_SuperLayers_FromStereo_ToAxial(vector<DCDCTrackCircle*>& locCDCTrackCircles, unsigned int locOuterSuperLayer, unsigned int locInnerSuperLayer);
      void Link_SuperLayers_FromStereo_ToStereo(vector<DCDCTrackCircle*>& locCDCTrackCircles, unsigned int locOuterSuperLayer, unsigned int locInnerSuperLayer);
      bool Check_IfShouldAttemptLink(const DCDCSuperLayerSeed* locSuperLayerSeed, bool locInnerSeedFlag);
      bool Attempt_SeedLink(DCDCSuperLayerSeed* locSuperLayerSeed1, DCDCSuperLayerSeed* locSuperLayerSeed2);
      bool Attempt_SeedLink(DCDCRingSeed& locRingSeed1, DCDCRingSeed& locRingSeed2, wire_direction_t locWireDirection1, wire_direction_t locWireDirection2);
      void Recycle_DCDCTrackCircle(DCDCTrackCircle* locCDCTrackCircle);

      // Continue DCDCTrackCircle Creation
      void Print_TrackCircles(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Print_TrackCircle(DCDCTrackCircle* locCDCTrackCircle);
      void Reject_DefiniteSpiralArms(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Drop_IncompleteGroups(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Fit_Circles(vector<DCDCTrackCircle*>& locCDCTrackCircles, bool locFitOnlyIfNullFitFlag, bool locAddStereoLayerIntersectionsFlag, bool locFitDuringLinkingFlag = false);
      DVector3 Find_IntersectionBetweenSuperLayers(const DCDCSuperLayerSeed* locInnerSuperLayerSeed, const DCDCSuperLayerSeed* locOuterSuperLayerSeed);

      // Filter Track Circles and Stereo Wires
      void Handle_StereoAndFilter(vector<DCDCTrackCircle*>& locCDCTrackCircles, bool locFinalPassFlag);
      void Truncate_TrackCircles(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Set_HitBitPattern_Axial(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Filter_TrackCircles_Axial(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Create_NewCDCSuperLayerSeeds(DCDCTrackCircle* locCDCTrackCircle);
      DCDCSuperLayerSeed* Create_NewStereoSuperLayerSeed(DCDCSuperLayerSeed* locCDCSuperLayerSeed, const DCDCTrackCircle* locCDCTrackCircle, map<DCDCTrkHit*, DCDCTrkHit*>& locProjectedStereoHitMap);
      bool Calc_StereoPosition(const DCDCWire *wire, const DHelicalFit* fit, DVector3 &pos, double &var_z, double& locPhiStereo, double d = 0.0);

      // Select Best Stereo Seeds and Calculate Theta & Z
      bool Select_CDCSuperLayerSeeds(DCDCTrackCircle* locCDCTrackCircle, bool locFinalPassFlag);
      void Select_ThetaZStereoHits(const DCDCTrackCircle* locCDCTrackCircle, int locInnerSeedSeriesIndex, int locOuterSeedSeriesIndex, bool locFinalPassFlag, vector<DCDCTrkHit*>& locComboHits);
      void Calc_StereoHitDeltaPhis(unsigned int locSuperLayer, vector<DCDCTrkHit*>& locHits, const DCDCTrackCircle* locCDCTrackCircle, vector<pair<DCDCTrkHit*, double> >& locDeltaPhis);
      double MinDeltaPhi(const vector<DCDCTrkHit*>& locInnerSeedHits, const vector<DCDCTrkHit*>& locOuterSeedHits);
      void Recycle_DCDCSuperLayerSeed(DCDCSuperLayerSeed* locCDCSuperLayerSeed);
      bool Find_ThetaZ(const DHelicalFit* locFit, const vector<DCDCTrkHit*>& locStereoHits, double& locTheta, double& locZ, double& locChiSqPerNDF);
      bool Find_ThetaZ_Regression(const DHelicalFit* locFit, const vector<DCDCTrkHit*>& locStereoHits, double& locTheta, double& locZ, double& locChiSqPerNDF);
      bool Find_Theta(const DHelicalFit* locFit, const vector<DCDCTrkHit*>& locStereoHits, double& locTheta, double& locThetaMin, double& locThetaMax, double& locChiSqPerNDF);
      bool Find_Z(const DHelicalFit* locFit, const vector<DCDCTrkHit*>& locStereoHits, double locThetaMin, double locThetaMax, double& locZ);
      void Set_HitBitPattern_All(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Filter_TrackCircles_Stereo(vector<DCDCTrackCircle*>& locCDCTrackCircles);

      // Finalize and Create DTrackCandidate Objects
      void Add_UnusedHits(vector<DCDCTrackCircle*>& locCDCTrackCircles);
      void Create_TrackCandidiate(DCDCTrackCircle* locCDCTrackCircle);
      bool Calc_PositionAndMomentum(DCDCTrackCircle* locCDCTrackCircle, DVector3 &pos, DVector3 &mom);

      int DEBUG_LEVEL;
      unsigned int MAX_ALLOWED_CDC_HITS;
      unsigned int MAX_NUM_RINGSEED_RINGS_SKIPABLE;
      double MAX_HIT_DIST; // cm
      double MAX_HIT_DIST2; // cm
      unsigned int MIN_SEED_HITS;
      unsigned int MAX_STRAWS_BETWEEN_LINK_SPIRAL_TURN;
      unsigned int MIN_STRAWS_POTENTIAL_SPIRAL_TURN;
      unsigned int MIN_STRAWS_DEFINITE_SPIRAL_TURN;
      unsigned int MIN_STRAWS_ADJACENT_TO_SPIRAL_TURN;
      unsigned int SEED_DENSITY_BIN_STRAW_WIDTH;
      unsigned int MAX_SEEDS_IN_STRAW_BIN;
      unsigned int MAX_ALLOWED_TRACK_CIRCLES; //bail if exceed this many
      bool ENABLE_DEAD_HV_BOARD_LINKING;
      unsigned int MAX_SUPERLAYER_NEW_TRACK; //don't allow new track seeds to start after this super layer //track seeds could start late if track is a decay product, or HV board is dead
      double MAX_COMMON_HIT_FRACTION;
      double MIN_CIRCLE_ASYMMETRY;
      double MAX_DRIFT_TIME;
      double MAX_SEED_TIME_DIFF;
      unsigned int MIN_PRUNED_STEREO_HITS;
      double MAX_UNUSED_HIT_LINK_ANGLE;

      size_t MAX_DCDCTrkHitPoolSize;
      size_t MAX_DCDCSuperLayerSeedPoolSize;
      size_t MAX_DCDCTrackCirclePoolSize;
      size_t MAX_HelicalFitPoolSize;

      double TARGET_Z;
      double VERTEX_Z_MIN;
      double VERTEX_Z_MAX;

      map<DCDCTrkHit*, unsigned int> dStereoHitNumUsedMap; //map of circle-projected (hit-z-group) hit to # DCDCSuperLayerSeeds it's in (when drops to 0, can recycle the hit memory)

      // Resource Pools for saving time and re-using memory
      vector<DCDCTrkHit*> dCDCTrkHitPool_Available;
      vector<DCDCTrkHit*> dCDCTrkHitPool_All;

      vector<DCDCSuperLayerSeed*> dCDCSuperLayerSeedPool_Available;
      vector<DCDCSuperLayerSeed*> dCDCSuperLayerSeedPool_All;

      vector<DCDCTrackCircle*> dCDCTrackCirclePool_Available;
      vector<DCDCTrackCircle*> dCDCTrackCirclePool_All;

      vector<DHelicalFit*> dHelicalFitPool_All;
      vector<DHelicalFit*> dHelicalFitPool_Available;

      vector<unsigned int> dNumStrawsPerRing; //index is ring index
      vector<unsigned int> superlayer_boundaries;

      unsigned int dNumCDCHits;
      vector<DCDCTrkHit*> cdctrkhits;
      vector<vector<DCDCTrkHit*> > cdchits_by_superlayer;
      vector<vector<DCDCSuperLayerSeed*> > dSuperLayerSeeds; //index 0 -> 6 is super layer 1 -> 7

      const DMagneticFieldMap* dMagneticField;
      double dFactorForSenseOfRotation;

      unsigned int dNumSeedDensityPhiBins;
      //dRejectedPhiRegions: due to hit density being too high in that region
      map<unsigned int, vector<pair<double, double> > > dRejectedPhiRegions; //key is super layer, pair is first/last phi (not necessarily min/max! (could pass thorugh phi = 0 barrier))
};

#endif // _DTrackCandidate_factory_CDC_