GlastFit.h Source File

00001 //------------------------------------------------------------------------------
00002 /* $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/reconstruction/GlastFit.h,v 1.1.1.1 1999/12/28 01:36:02 burnett Exp $ */
00003 // File: G
00004 // Contents ----------------------------------------------------------------
00005 //
00006 //     GlastFit
00007 //
00008 //             Is a Kalman Track Follower class
00009 //
00010 //             It uses a Ray direction as starting seed
00011 //             and picks Si-Clusters inside a n-sigma region of the projected track into the Si Plane
00012 //
00013 //
00014 //                              B. Atwood
00015 //                              B. Atwood, JA Hernando    Santa Cruz 2/26/99
00016 //------------------------------------------------------------------------------
00017 
00018 #ifndef __GLASTFIT_H
00019 #define __GLASTFIT_H 1
00020 
00021 //#include "TrackerRecon.h"
00022 #include "KalFit.h"
00023 #include "geometry/Ray.h"
00024 
00025 class Ray;
00026 class SiClusters;
00027 //----------------------------------------
00028 //   GlastFit
00029 //----------------------------------------
00030 //################################
00031 class GlastFit : public KalTrack
00032 //################################
00033 {
00034 public:
00035     
00036     GlastFit(enum SiData::Axis a, SiClusters*, float cut, float energy,
00037         int ist, const Ray& testRay);
00038     virtual ~GlastFit() {}
00039 
00040     // Original Glast Fit functions - Compatibility
00041     //----------------------------------------------
00042     float truncateFitAfter (int layer);
00043     void flagAllHits ();
00044     void unFlagAllHits ();
00045     void flagHit (int idata);
00046     void unFlagHit (int idata);
00047     
00048     int seedHits();
00049     double seedChiSq(int nhits);
00050 
00051     // Acces functions
00052     //-------------------
00053     SiData::Axis getAxis() { return m_axis;}
00054     int gaps ()  {  return m_gaps;}
00055     int firstGaps () { return m_istGaps;}
00056     int noise() { return m_total_noise;}
00057     int startNoise() { return m_ist_noise;}
00058     int firstLayer() { return m_istLayer;}
00059     int lastLayer() { return m_lstLayer;}
00060     float energy () { return m_iniEnergy;}
00061     float Q() { return m_quality;}
00062 
00063 protected:
00064     
00065     double getZklayer(enum SiData::Axis axis, int klayer);
00066     bool crack(const KalPlane&);
00067     Point doVertex(const GlastFit&);
00068     
00069 private:
00070     
00071     // Initial information
00072     //---------------------
00073     SiData::Axis m_axis;
00074     double m_sigmaCut;
00075     float m_iniEnergy;
00076     int m_iniLayer;
00077     SiClusters* _mdata;
00078 
00079     // Contabability
00080     //----------------
00081     int m_gaps; 
00082     int m_istGaps;
00083     int m_istLayer;
00084     int m_lstLayer;
00085     int m_total_noise;
00086     int m_ist_noise;
00087     int m_lstGaps;
00088 
00089         double m_quality;
00090   
00091 };
00092 
00093 //---------------------------------------------------------------
00094 //
00095 //    GF - Tracking classes
00096 //
00097 //                     B.Atwood, JA Hernando
00098 //                     Santa Cruz, 03/30/99
00099 //---------------------------------------------------------------
00100 
00101 //############################################
00102 class GFdata
00103 //############################################
00104 {
00105 public:
00106 
00107     GFdata() : m_RCenergy(0.), m_quality(0.), m_firstLayer(0), 
00108                m_nhits(0), m_itower(0)
00109     {ini();}
00110                 
00111     GFdata getGFdata() {return *this;}
00112     
00113         // access
00114     inline Point vertex() const {return m_vertex;}
00115     inline Vector direction() const {return m_direction.unit();}
00116     inline double RCenergy() const {return m_RCenergy;}
00117     inline double Q() const {return m_quality;}
00118     inline int firstLayer() const {return m_firstLayer;}
00119     inline int nhits() const {return m_nhits;}
00120     inline int tower() const {return m_itower;}
00121 
00122     Ray ray() const {return Ray(m_vertex,m_direction);} 
00123 
00124         // operations
00125         void ini();
00126         bool empty() const;
00127         void writeOut(int level) const;
00128 
00129         // utilities
00130     static Point doVertex(const Ray&, const  Ray& );
00131         static Vector doDirection(const Vector& xdir, const Vector& ydir);
00132     static bool neighbourTowers(int itower, int jtower);
00133 
00134 public:
00135 
00136         enum particle { ELECTRON, NOELECTRON };
00137         enum type {TRACK, PARTICLE, PAIR, GAMMA};
00138 
00139 protected:
00140 
00141         // Output Data
00142     Point m_vertex;
00143     Vector m_direction;
00144     double m_RCenergy;
00145     double m_quality;
00146     int m_firstLayer;
00147         int m_nhits;
00148     int m_itower;
00149         
00150 };
00151  
00152 //############################################
00153 class GFbase: public GFdata
00154 //############################################
00155 {
00156         
00157 
00158 public:
00159 
00160         // access
00161         inline double sigmaCut() const {return m_sigmaCut;}
00162         inline bool alive() const {return m_alive;}
00163         inline Point inputVertex() const {return m_inVertex;}
00164         inline Vector inputDirection() const {return m_inDirection;}
00165         inline int inputLayer() const {return m_iniLayer;}
00166         inline double inputEnergy() const {return m_iniEnergy;}
00167 
00168         inline Ray inputRay() const {return Ray(m_inVertex,m_inDirection);}
00169 
00170         // operation
00171         virtual void flagAllHits() = 0;
00172         virtual void unFlagAllHits() = 0;
00173 
00174         virtual bool empty() const = 0;
00175         virtual bool accept() const = 0;
00176         virtual void clear() = 0;
00177         virtual void writeOut(int level) const = 0; 
00178 
00179 
00180 public:
00181         
00182     enum StatusHit {EMPTY, FOUND, CRACK};
00183     enum StatusPair {TOGETHER, SPLIT, ONE, DONE};
00184 
00185 protected:
00186 
00187         // constructor
00188     GFbase(double sigmaCut, double ene, int ist, const Ray& testRay);
00189     virtual ~GFbase() {}
00190 
00191         // operations
00192         virtual void ini() = 0;
00193         virtual void doit(); // do a full pattern Regonition
00194         virtual void step(int kplane) = 0; // One Step in the Pattern Recognition
00195         virtual void anastep(int kplane) = 0; // make some analysis after the step is done
00196         virtual void fit() = 0; // fit the GF object - compute the GFdata
00197         virtual bool end() const = 0; // end of the Pattern Recognition?
00198         virtual void kill() = 0;        
00199         virtual void setAlive() = 0;
00200 
00201         virtual void contability(int kplane) = 0; // contability of the anaStep
00202         virtual void loadGFdata() = 0;   // load the GF data
00203 
00204 protected:
00205     
00206         // Control 
00207     double m_sigmaCut;
00208         bool m_alive;
00209 
00210         // Input Data
00211         Point m_inVertex;
00212         Vector m_inDirection;
00213     double m_iniEnergy;
00214     int m_iniLayer;
00215 
00216 };
00217 
00218 //##############################################
00219 class GFsegment: public KalTrack
00220 //##############################################
00221 {
00222 protected:
00223 
00224         friend class GFgamma;
00225         friend class GFparticle;
00226         friend class GFpair;
00227         friend class GFtrack;
00228 
00229         // set
00230         GFsegment(const GFtrack* _GFtrack);
00231 
00232         // access
00233         int indexhit() const {return m_indexhit;}
00234         GFbase::StatusHit status() const {return m_statusHit;}
00235         KalPlane getKPlane() const;
00236         double chiGFSq() const;
00237 
00238 
00239         // operations
00240         void best(int kplane);
00241         void next(int kplane);
00242         void previous(int kplane); 
00243         void clear();
00244         bool accept() const;
00245 
00246         void flagUsedHits(int kplane);
00247         void unFlagAllHits();
00248 
00249 private:
00250             
00251         // Project to the next Plane
00252         KalPlane followingKPlane(int kplane) const;
00253         KalPlane getKPlane(int kplane) const;
00254         void doit(KalPlane& oriKplane, int jplane, KalHit::TYPE type = KalHit::FIT);
00255     KalPlane projectedKPlane(KalPlane previous, int klayer, 
00256                 KalHit::TYPE type = KalHit::FIT) const;
00257         GFbase::StatusHit nextKPlane(const KalPlane& previous, int kplane, 
00258                 KalPlane& next, KalHit::TYPE typ = KalHit::FIT) const; // returns next n
00259     
00260         // Selecting the Hit
00261     double sigmaFoundHit(const KalPlane& previous, const KalPlane& next, int& indexhit, double& radius) const; // returns also indexhit and radius
00262         void incorporateFoundHit(KalPlane& next, int indexhit) const; // modifies next
00263     bool foundHit(int& indexhit, double& inerRadius, double outRadius,
00264                         const Point& CenterX, const Point& nearHit, double slope) const;
00265 
00266     // Utilities
00267     double getZklayer(enum SiData::Axis axis, int klayer) const;
00268     bool crack(const KalPlane&) const;
00269 
00270 private:
00271 
00272         SiData::Axis m_axis;
00273 
00274         KalPlane m_nextKplane;
00275         int m_indexhit;
00276         GFbase::StatusHit m_statusHit;
00277 
00278         const GFtrack* _mGFtrack;
00279 
00280 };
00281 
00282 //##############################################
00283 class GFtrack: public GFbase, public KalTrack
00284 //##############################################
00285 {    
00286 public:
00287     
00288     GFtrack(enum SiData::Axis axis, double sigmaCut,
00289                double energy, 
00290                int ist, 
00291                const Ray& testRay, bool doit = true);
00292         ~GFtrack() {
00293                 delete _mGFsegment;
00294         }
00295         
00296         //-- from GFdata 
00297         void flagAllHits();
00298         void unFlagAllHits();
00299         
00300         bool empty() const;
00301         bool accept() const;
00302         void clear();
00303         void writeOut(int level) const; 
00304         //--
00305 
00306     // acces 
00307     SiData::Axis getAxis() const {return m_axis;}
00308     int numGaps ()  const {return m_gaps;}
00309     int numFirstGaps () const {return m_istGaps;}
00310     int numNoise() const {return m_noisyHits;}
00311         int numFirstNoise() const {return m_istNoisyHits;}
00312     int lastLayer() const { return m_lstLayer;}
00313 
00314         // operations
00315         bool veto(int& indexhit, double& sigma) const;
00316         double Qbest() const {return m_qbest;};
00317         
00318         double computeQuality() const;
00319 
00320 protected:
00321 
00322         friend class GFsegment;
00323 
00324         friend class GFparticle;
00325         friend class GFpair;
00326         friend class GFgamma;
00327 
00328         //-- from GFbase
00329         void ini();
00330 //      void doit(); // do a full pattern Regonition
00331         void step(int kplane); // One Step in the Pattern Recognition
00332         void anastep(int kplane); // make some analysis after the step is done
00333         void fit(); // fit the GF object - compute the GFdata
00334         bool end() const; // end of the Pattern Recognition?
00335         void kill();    
00336         void setAlive();
00337 
00338         void contability(int kplane); // contability of the anaStep
00339         void loadGFdata();   // after the fit load the data
00340         // ---
00341 
00342         // creation
00343         void setIniEnergy(double ene);
00344         void setStatus(StatusHit status) {m_status = status;}
00345 
00346         // access to the Step Plane 
00347         StatusHit status() const {return m_status;}
00348         KalPlane firstKPlane() const;
00349         KalPlane lastKPlane() const;
00350         KalPlane previousKPlane() const;
00351     KalPlane originalKPlane() const;
00352 
00353         // operations
00354         void removeStep(int kplane = -1);
00355         double doQbest();
00356         void associateOrthStep(const GFtrack* _OrhGFtrack, KalHit::TYPE type = KalHit::FIT);
00357         void associateOrthGFtrack(const GFtrack* _OrhGFtrack, bool fix = false, 
00358                 KalHit::TYPE type = KalHit::FIT);
00359 
00360 private:
00361     
00362         // Parasitus 
00363         GFsegment* _mGFsegment;
00364 
00365     // Input Data
00366     SiData::Axis m_axis;
00367 
00368     // Status
00369         StatusHit m_status;
00370         int m_lstGaps;
00371 
00372         // Output Data
00373         double m_qbest;
00374 
00375         // contability
00376     int m_gaps; 
00377     int m_istGaps;
00378     int m_lstLayer;
00379     int m_noisyHits;
00380     int m_istNoisyHits;
00381 };
00382 
00383 //##############################################
00384 class GFparticle: public GFbase
00385 //##############################################
00386 {    
00387 public:
00388     
00389     GFparticle(double sigmaCut, double energy, int ist, 
00390                const Ray& testRay, bool doit = true);
00391         ~GFparticle() {
00392                 delete _mXGFtrack;
00393                 delete _mYGFtrack;
00394         }
00395         
00396         //-- from GFdata 
00397         void flagAllHits();
00398         void unFlagAllHits();
00399         
00400         bool empty() const;
00401         bool accept() const;
00402         void clear();
00403         void writeOut(int level) const; 
00404         //--
00405 
00406     // access 
00407         const GFtrack* getXGFtrack() const {return _mXGFtrack;}
00408         const GFtrack* getYGFtrack() const {return _mYGFtrack;}
00409 
00410     int numGaps ()  const {return m_gaps;}
00411     int numFirstGaps () const {return m_istGaps;}
00412     int numNoise() const {return m_noisyHits;}
00413         int numFirstNoise() const {return m_istNoisyHits;}
00414     int lastLayer() const { return m_lstLayer;}
00415 
00416         // operations
00417         bool veto(int& indexhit, double& sigma) const;
00418         double Qbest() const {return m_qbest;};
00419 
00420 protected:
00421 
00422         friend class GFsegment;
00423 
00424         // friend class GFparticle;
00425 
00426         friend class GFpair;
00427         friend class GFgamma;
00428 
00429         //-- from GFbase
00430         void ini();
00431 //      void doit(); // do a full pattern Regonition
00432         void step(int kplane); // One Step in the Pattern Recognition
00433         void anastep(int kplane); // make some analysis after the step is done
00434         void fit(); // fit the GF object - compute the GFdata
00435         bool end() const; // end of the Pattern Recognition?
00436         void kill();    
00437         void setAlive();
00438 
00439         void contability(int kplane); // contability of the anaStep
00440         void loadGFdata();   // after the fit load the data
00441         // ---
00442 
00443         // creation
00444         void setIniEnergy(double ene);
00445         void setStatus(StatusHit status) {m_status = status;}
00446 
00447         // access
00448         StatusHit status() const {return m_status;}     
00449 
00450         // association operations
00451         void associateStatus();
00452         void associateStep();
00453         void associateAnaStep();
00454         void associateFit();
00455 
00456         // operations
00457         double doQbest();
00458 
00459         // utilities
00460         static bool sameTower(const GFtrack* _GFtrk1,const GFtrack* _GFtrk2);
00461         static bool removeWorseStep(GFtrack* _GFtrkX, GFtrack* _GFtrkY);
00462 
00463 private:
00464     
00465     // Status
00466         bool m_associate;   // Do a 3D track
00467         bool m_conflictPattern; // Topology conflict 
00468         StatusHit m_status;
00469 
00470         // Output Data
00471         double m_qbest;
00472 
00473         // contability
00474     int m_gaps; 
00475     int m_istGaps;
00476     int m_lstLayer;
00477     int m_noisyHits;
00478     int m_istNoisyHits;
00479 
00480         // tracks
00481         GFtrack* _mXGFtrack;
00482         GFtrack* _mYGFtrack;
00483 };
00484 
00485 //#############################################################
00486 class GFpair : public GFbase
00487 //#############################################################
00488 {
00489 public:   
00490 
00491         // construction
00492     GFpair(double xene, enum SiData::Axis axis, double sigmaCut,
00493                 double energy,int ist, const Ray& testRay,
00494                 bool doit = true);
00495     ~GFpair() { 
00496                 delete _mGFbest;
00497                 delete _mGFpair;
00498     }
00499 
00500         //-- from GFdata 
00501         void flagAllHits();
00502         void unFlagAllHits();
00503         
00504         bool empty() const;
00505         bool accept() const;
00506         void clear();
00507         void writeOut(int level) const; 
00508         //--
00509 
00510     // access 
00511     GFtrack* getBest() const {return _mGFbest;}
00512     GFtrack* getPair() const {return _mGFpair;}
00513 
00514         double weightSlope() const {return m_weightBest;}
00515         double errorSlope() const {return m_errorSlope;}
00516 
00517         int numTogether() const {return m_together;}
00518         int numSplit() const {return m_split;}
00519         int numOne() const {return m_one;}
00520         int numSharedHits() const {return m_shared;}
00521         int numEmpty() const {return m_empty;}
00522         
00523 protected:
00524 
00525         //-- from GFbase
00526         void ini();
00527 //      void doit(); // do a full pattern Regonition
00528         void step(int kplane); // One Step in the Pattern Recognition
00529         void anastep(int kplane); // make some analysis after the step is done
00530         void fit(); // fit the GF object - compute the GFdata
00531         bool end() const; // end of the Pattern Recognition?
00532         void kill();    
00533         void setAlive();
00534 
00535         void contability(int kplane); // contability of the anaStep
00536         void loadGFdata();   // after the fit load the data
00537         // ---
00538 
00539         // Set
00540         void setIniEnergy();
00541         void setDecideBest(bool decideBest) {m_decideBest = decideBest;}
00542         void setStatus(StatusPair newStatus);
00543 
00544         // Status
00545         StatusPair status() const {return m_status;}
00546     void newStatus(int klayer);
00547         bool forceSplit(int klayer) const; 
00548 
00549         // Stepping
00550     void stepTogether(int kplane);
00551     void stepSplit(int kplane);
00552         void selfishStepSplit(int kplane);
00553 
00554         // Selecting the Best Track
00555         void decideBest();
00556         void swap();
00557 
00558         // GFdata - Utilities
00559         Vector doDirection(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2,
00560                 double& weight1, double& errorSlope);
00561         Vector doDirection(double& weight1);
00562         Vector doDirectionXene(double xene, double& weight1);
00563         double doEnergy(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2);
00564 
00565         // Stepping - Utilities
00566         bool allowedShareHit(const GFtrack* _GFtrack) const;
00567         void removeWorseStep(GFtrack* _GFtrk1, GFtrack* _GFtrk2); 
00568         void resizeSharedHits();
00569 
00570 private:
00571     
00572     friend class GFgamma;
00573     // in data
00574     double m_xEne;
00575     SiData::Axis m_axis;
00576 
00577     StatusPair m_status;
00578     bool m_decideBest;
00579 
00580         // Slope determination
00581         double m_weightBest;
00582         double m_errorSlope;
00583 
00584     // contability
00585     int m_together;
00586     int m_split;
00587     int m_one;
00588     int m_shared;
00589     int m_empty;
00590   
00591 protected:
00592 
00593         // output data
00594     GFtrack* _mGFbest;
00595     GFtrack* _mGFpair;
00596 
00597     GFtrack* _mGFalive;
00598 };
00599 
00600 //############################################
00601 class GFgamma : public GFbase
00602 //############################################
00603 {
00604 public:
00605 
00606     GFgamma(double xene,
00607                 double sigmaCut,
00608                 double energy,  
00609                 int ist, 
00610                 const Ray& testRay);
00611     ~GFgamma() {
00612                 delete _mXpair;
00613                 delete _mYpair;
00614     }
00615 
00616         //-- from GFdata 
00617         void flagAllHits();
00618         void unFlagAllHits();
00619         
00620         bool empty() const;
00621         bool accept() const;
00622         void clear();
00623         void writeOut(int level) const; 
00624         //--
00625 
00626         // Access
00627         bool conflictPattern() const {return m_conflictPattern;}
00628         bool fix() const {return m_fixTopology;}
00629     GFpair* getXpair() const {return _mXpair;}
00630     GFpair* getYpair() const {return _mYpair;}
00631         GFtrack* getBest(SiData::Axis axis) const {
00632                 if (axis == SiData::X) return _mXpair->getBest();
00633                 else return _mYpair->getBest();
00634         }
00635         GFtrack* getPair(SiData::Axis axis) const {
00636                 if (axis == SiData::X) return _mXpair->getPair();
00637                 else return _mYpair->getPair();
00638         }
00639         Point getFirstHit() const;
00640 
00641         int numTogether() const {return m_together;}
00642         int numSplit() const {return m_split;}
00643         int numOne() const {return m_one;}
00644 
00645         // operations
00646         bool veto() const;
00647         double Qbest();
00648     static bool accept(const GFdata&, const GFdata&);
00649 
00650 protected:
00651 
00652         //-- from GFbase
00653         void ini();
00654 //      void doit(); // do a full pattern Regonition
00655         void step(int kplane); // One Step in the Pattern Recognition
00656         void anastep(int kplane); // make some analysis after the step is done
00657         void fit(); // fit the GF object - compute the GFdata
00658         bool end() const; // end of the Pattern Recognition?
00659         void kill();    
00660         void setAlive();
00661 
00662         void contability(int kplane); // contability of the anaStep
00663         void loadGFdata();   // after the fit load the data
00664         // ---
00665 
00666         // set
00667         void construct();
00668         void setDecideBest(bool decideBest);
00669         
00670         // Stepping
00671         StatusPair newStatus();
00672         void connectStep();
00673         void associateStep();
00674         void topologyStep();
00675         void associateStatus(StatusPair status);
00676 
00677         void associateAnaStep();
00678         void associateAnaStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2);
00679 
00680         // Best tracks
00681         void decideBest();
00682         void associateFit();
00683 
00684         // Utilities
00685         static bool crossingTowers(const GFtrack* _Xtrk1, const GFtrack* _Ytrk1,
00686                 const GFtrack* _Xtrk2, const GFtrack* _Ytrk2);
00687 
00688 private:
00689 
00690     double m_xEne;
00691         // status
00692         // The Status is check by errorControl(); - some combination are not allowed!
00693         
00694         bool m_connect; // make a 3D gamma or a 2D gamma
00695         bool m_associate;  // association of the 2D pairs (need patternSwap to be defined!)
00696         bool m_patternSwap; // false (best-best), true (best-pair) the two different associations
00697         bool m_fixTopology; // the topology is fixed
00698         bool m_decideBest; // the decition of the best track
00699         
00700         bool m_conflictPattern; // true is a pattern Conflic is detected and can not be solved
00701         bool m_swapDone;       // Swapping completed
00702 
00703         // Status
00704         StatusPair m_status;
00705 
00706         // contability
00707     int m_together;
00708         int m_split;
00709         int m_one;
00710 
00711     GFpair* _mXpair;
00712     GFpair* _mYpair;
00713 };
00714 
00715 //############################################
00716 class GFtutor
00717 //############################################
00718 {
00719 public:
00720 
00721         static SiClusters* _DATA;
00722 
00723         static bool CONTROL_connectGFpair;
00724         static bool CUT_veto;
00725 
00726 public:
00727 
00728     static int  okClusterSize(SiData::Axis axis, int indexhit, double slope);   
00729     static bool neighbourTowers(int itower, int jtower);
00730 }; 
00731 
00732 #endif