FastTracklet.h

Go to the documentation of this file.

/*

Supporting structure of fast tracking, include SignedHit and Tracklet

Author: Kun Liu, liuk@fnal.gov
Created: 06-09-2013

*/

#ifndef _FASTTRACKLET_H
#define _FASTTRACKLET_H

#include <GlobalConsts.h>

#include <list>
#include <vector>

#include <TVector3.h>

#include <geom_svc/GeomSvc.h>
#include <phool/PHObject.h>

#include "SRawEvent.h"
#include "SRecEvent.h"

class SignedHit : public PHObject
{
public:
    SignedHit();
    explicit SignedHit(int detectorID);
    SignedHit(Hit hit_input, int sign_input);

    //PHObject virtual overloads
    void identify(std::ostream& os = std::cout) const;
    void Reset() { hit = Hit(); sign = 0; }
    int  isValid() const { return hit.index > 0; }
    SignedHit* Clone() const { return (new SignedHit(hit, sign)); }

    //comparision operators for sorting
    bool operator<(const SignedHit elem) const { return hit.detectorID < elem.hit.detectorID; }
    bool operator==(const SignedHit elem) const { return hit.index == elem.hit.index; }

    //Get the real hit position
    double pos() { return hit.pos + sign*hit.driftDistance; }
    double pos(int sign_input) { return hit.pos + sign_input*hit.driftDistance; }

    //Data members
    Hit hit;
    int sign;

    ClassDef(SignedHit, 1)
};

class PropSegment : public PHObject
{
public:
    PropSegment();

    //PHObject virtual overloads
    void identify(std::ostream& os = std::cout) const { print(os); };
    void Reset() { init(); }
    int  isValid() const;
    PropSegment* Clone() const { return (new PropSegment(*this)); }

    //init -- temporary, only used for tests
    void init();

    //Quality cut
    //bool isValid();

    //Debugging output
    void print(std::ostream& os = std::cout) const;

#ifndef __CINT__
    //Get expected position at a given Z
    double getExpPosition(double z) const { return a*z + b; }

    //Get the closest approach to a given space position/proptube wire
    double getClosestApproach(double z, double pos);

    //Get reference pos at first two planes
    double getPosRef(double default_val = -9999.);

    //Number of hits
    int getNHits() const;

    //Number of planes
    int getNPlanes() const;

    //Fit the segment -- naive linear fit
    void fit();   // external call
    void fit_2hits();
    void fit_34hits();

    //linear chisq fitter
    void linearFit_simple();
    void linearFit_iterative();

    //resolve left/right
    void resolveLR();
    void resolveLR(int setting);
#endif

    //track slope the interception
    double a;
    double b;
    double err_a;
    double err_b;

    //chisq of the segment
    double chisq;

    //Auxilary hodoscope hit list, cannot be possibly more than 4 hodoscope hits
    int nHodoHits;
    Hit hodoHits[4];

    //Hit list -- only 4 hits at most
    SignedHit hits[4];

    ClassDef(PropSegment, 4)
};

class Tracklet : public PHObject
{
public:
    Tracklet();

    //PHObject virtual overloads
    void identify(std::ostream& os = std::cout) const { os << "Tracklet @sID=" << stationID << std::endl;}
    void Reset() { *this = Tracklet(); }
    int  isValid() const;
    Tracklet* Clone() const { return (new Tracklet(*this)); }

    //Basic quality cut
    //bool isValid();

    //Debuggin output
    void print(std::ostream& os = std::cout);

#ifndef __CINT__
    //Sort hit list
    void sortHits() { hits.sort(); }

    //Update/get number of real hits
    void updateNHits();
    int getNHits() const { return nXHits + nUHits + nVHits; }

    //Number of all hits (even excluded)
    int getNAllHits() { return hits.size(); }

    //Get the probabilities
    double getProb() const;

    //Get the momentum probabilities
    double getMomProb() const;

    //Get the chi square
    double getChisq() const { return chisq; }

    //Get x and y positions at a given z
    double getExpPositionX(double z) const;
    double getExpPosErrorX(double z) const;
    double getExpPositionY(double z) const;
    double getExpPosErrorY(double z) const;
    double getExpPositionW(int detectorID);

    //Get momentum upstream/downstream
    TVector3 getMomentumSt1();
    TVector3 getMomentumSt3();
    TVector3 getExpMomentum(double z);

    //Get the i-th signed hit
    SignedHit getSignedHit(int index);

    //Kernal function to calculate chi square for minimizer
    double Eval(const double* par);
    double calcChisq();

    //Add dummy hits
    void addDummyHits();

    //Momentum estimation using back partial
    double getMomentum() const;

    //Decide charge by KMag bending direction
    //int getCharge() const { return x0*KMAGSTR > tx ? 1 : -1; }
    int getCharge() const;

    //Get the slope and intersection in station 1
    void getXZInfoInSt1(double& tx_st1, double& x0_st1);
    void getXZErrorInSt1(double& err_tx_st1, double& err_x0_st1);

    //For sorting tracklet list
    bool operator<(const Tracklet& elem) const;

    //For reducing similar tracklets
    bool similarity(const Tracklet& elem) const;

    //Merge the hit list from two tracklets
    Tracklet merge(Tracklet& elem);

    //For adding two tracklets together to form a back partial track
    Tracklet operator+(const Tracklet& elem) const;

    //For adding two tracklets together to form a global track
    Tracklet operator*(const Tracklet& elem) const;

    //Convert to a SRecTrack
    SRecTrack getSRecTrack(bool hyptest = true);
#endif

    //Station ID, ranging from 1 to nStation, nStation-1 means back partial track, nStation means global track
    int stationID;

    //Number of hits
    int nXHits;
    int nUHits;
    int nVHits;

    //Chi square
    double chisq;

    //chisq at vertex
    double chisq_vtx;

    //List of signed hits
    std::list<SignedHit> hits;

    //Corresponding prop. tube segments
    PropSegment seg_x;
    PropSegment seg_y;

    //Slope, intersection, momentum and their errors
    double tx;
    double ty;
    double x0;
    double y0;
    double invP;

    double err_tx;
    double err_ty;
    double err_x0;
    double err_y0;
    double err_invP;

    //Residuals of all pos
    double residual[nChamberPlanes];

    ClassDef(Tracklet, 4)
};

class TrackletVector : public PHObject
{
public:
    TrackletVector();
    virtual ~TrackletVector();

    void identify(std::ostream& os = std::cout) const;
    void Reset();
    int  isValid() const { return 1; };
    TrackletVector* Clone() const { return (new TrackletVector(*this)); }

    bool empty() const { return trackletVec.empty(); }
    size_t size() const { return trackletVec.size(); }
    void clear() { Reset(); }

    const Tracklet* at(const size_t index) const;
    Tracklet* at(const size_t index);
    void push_back(const Tracklet* tracklet);
    size_t erase(const size_t index);

    std::vector<Tracklet*>::const_iterator begin() const { return trackletVec.begin(); }
    std::vector<Tracklet*>::const_iterator end()   const { return trackletVec.end(); }  

    std::vector<Tracklet*>::iterator begin() { return trackletVec.begin(); }
    std::vector<Tracklet*>::iterator end()   { return trackletVec.end(); }

private:
    std::vector<Tracklet*> trackletVec;

    ClassDef(TrackletVector, 1)
};


#endif