SiliconTracking_MarlinTrk.h

#ifndef SILICONTRACKING_MarlinTrk_H
#define SILICONTRACKING_MarlinTrk_H 1

#include "marlin/Processor.h"
#include <marlin/Global.h>
#include "lcio.h"
#include <string>
#include <vector>
#include <cmath>
#include <IMPL/TrackImpl.h>
#include "ClusterExtended.h"
#include "TrackExtended.h"
#include "TrackerHitExtended.h"
#include "HelixClass.h"

#include "MarlinTrk/IMarlinTrack.h"

#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>

namespace DiagnosticsHistograms {
  class Histograms ;
}

namespace DiagnosticsHistograms2D {
  class Histograms ;
}

using namespace lcio ;
using namespace marlin ;

namespace  gear{
  class GearMgr ;
}

namespace MarlinTrk {
  class HelixFit;
  class IMarlinTrkSystem ;
}

namespace UTIL{
  class LCRelationNavigator ;
}


class SiliconTracking_MarlinTrk : public Processor {
  
public:
  
  virtual Processor*  newProcessor() { return new SiliconTracking_MarlinTrk ; }
  
  
  SiliconTracking_MarlinTrk() ;
  
  virtual void init() ;
  
  virtual void processRunHeader( LCRunHeader* run ) ;
  
  virtual void processEvent( LCEvent * evt ) ; 
  
  
  virtual void check( LCEvent * evt ) ; 
  
  
  virtual void end() ;
  
  
protected:
  
  int _nRun ;
  int _nEvt ;
  EVENT::LCEvent* _current_event;
  
  int _nDivisionsInPhi;
  int _nDivisionsInTheta;
  int _nLayers;
  
  MarlinTrk::HelixFit* _fastfitter;
  
  MarlinTrk::IMarlinTrkSystem* _trksystem ;
  bool _runMarlinTrkDiagnostics;
  std::string _MarlinTrkDiagnosticsName;
  
  bool _MSOn, _ElossOn, _SmoothOn ;
  
  float _initialTrackError_d0;
  float _initialTrackError_phi0;
  float _initialTrackError_omega;
  float _initialTrackError_z0;
  float _initialTrackError_tanL;
  
  double _maxChi2PerHit;  
  
  bool  _UseEventDisplay;
  int _detector_model_for_drawing;
  std::vector<int> _colours;  
  float     _helix_max_r;
  
  void drawEvent();
  
  
  // histogram member variables
  
  bool  _createDiagnosticsHistograms;
  DiagnosticsHistograms::Histograms* _histos ;

  
  int _ntriplets, _ntriplets_good, _ntriplets_2MCP, _ntriplets_3MCP, _ntriplets_1MCP_Bad, _ntriplets_bad;
  
  
  LCCollection* GetCollection(  LCEvent * evt, std::string colName ) ;
  
  LCRelationNavigator* GetRelations( LCEvent * evt, std::string RelName ) ;
  
  std::string  _colNameMCParticles;
  float _MCpThreshold ;
  
  
  struct compare_TrackExtended{
    // n.b.: a and b should be TrackExtended const *, but the getters are not const :-(
    bool operator()(TrackExtended *a, TrackExtended *b) const {
      if ( a == b ) return false;
      return (a->getChi2()/a->getNDF() < b->getChi2()/b->getNDF() );
    }
  };
  
  
  std::string _VTXHitCollection;
  std::string _FTDPixelHitCollection;
  std::string _FTDSpacePointCollection;
  std::string _SITHitCollection;
  std::string _siTrkCollection;
  
  std::vector< LCCollection* > _colTrackerHits;
  std::map< LCCollection*, std::string > _colNamesTrackerHits;
  
  std::vector<TrackerHitExtendedVec> _sectors;
  std::vector<TrackerHitExtendedVec> _sectorsFTD;
  
  class TracksWithNHitsContainer {
  public:
    void clear();
    
    inline void resize(size_t maxHits) {
      _tracksNHits.resize(maxHits-2);
      _maxIndex=(maxHits-3);
    }
    
    // Sort all track vectors according to chi2/nDof
    //      void sort();
    
    inline TrackExtendedVec & getTracksWithNHitsVec( size_t nHits ) {
      //return _tracksNHits[ std::min(nHits-3, _maxIndex) ];
      // for debugging: with boundary check
      return _tracksNHits.at(std::min(nHits-3, _maxIndex));
    }
    
  protected:
    std::vector< TrackExtendedVec > _tracksNHits;
    size_t _maxIndex; 
  };
  
  TracksWithNHitsContainer _tracksWithNHitsContainer;
  
  int InitialiseVTX(LCEvent * evt);
  int InitialiseFTD(LCEvent * evt);
  void ProcessOneSector(int iSectorPhi, int iSectorTheta);
  void CleanUp();
  TrackExtended * TestTriplet(TrackerHitExtended * outerHit, 
                              TrackerHitExtended * middleHit,
                              TrackerHitExtended * innerHit,
                              HelixClass & helix);
  
  int BuildTrack(TrackerHitExtended * outerHit, 
                 TrackerHitExtended * middleHit,
                 TrackerHitExtended * innerHit,
                 HelixClass & helix, 
                 int innerlayer,
                 int iPhiLow, int iPhiUp,
                 int iTheta, int iThetaUp,
                 TrackExtended * trackAR);
  
  void Sorting( TrackExtendedVec & trackVec);
  void CreateTrack(TrackExtended * trackAR );
  void AttachRemainingVTXHitsSlow();
  void AttachRemainingFTDHitsSlow();
  void AttachRemainingVTXHitsFast();
  void AttachRemainingFTDHitsFast();
  void TrackingInFTD();
  int BuildTrackFTD(TrackExtended* trackAR, int* nLR, int iS);
  int AttachHitToTrack(TrackExtended * trackAR, TrackerHitExtended * hit, int iopt);
  
  void FinalRefit(LCCollectionVec* trk_col, LCCollectionVec* rel_col);
  
  float _bField;
  float _chi2WRPhiTriplet;
  float _chi2WRPhiQuartet;
  float _chi2WRPhiSeptet;
  float _chi2WZTriplet;
  float _chi2WZQuartet;
  float _chi2WZSeptet;
  float _minDistCutAttach;
  int _minimalLayerToAttach;
  
  // two pi is not a constant in cmath. Calculate it, once!
  static const double TWOPI;
  
  double _dPhi;
  double _dTheta;
  double _dPhiFTD;
  

  
  std::vector<int> _Combinations;
  std::vector<int> _CombinationsFTD;
  
  float _resolutionRPhiVTX;
  float _resolutionZVTX;
  
  float _resolutionRPhiFTD;
  float _resolutionZFTD;
  
  float _resolutionRPhiSIT;
  float _resolutionZSIT;
  
  float _phiCutForMerging;
  float _tanlambdaCutForMerging;
  float _angleCutForMerging;
  
  int _print;
  int _checkForDelta;
  float _minDistToDelta;
  
  float _distRPhi;
  float _distZ;
  float _chi2FitCut;
  
  
  TrackExtendedVec _trackImplVec;
  
  
  float _cutOnD0, _cutOnZ0, _cutOnOmega, _cutOnPt;
  
  int _minimalHits;
  int _nHitsChi2;
  int _attachFast;
  
  int _max_hits_per_sector;
  
  int _nTotalVTXHits,_nTotalFTDHits,_nTotalSITHits;
  int _useSIT;

  std::string _trkSystemName ;
  
  //  int _createMap;
  
  UTIL::BitField64* _encoder;
  int getDetectorID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::subdet]; }
  int getSideID(TrackerHit* hit)     { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::side]; };
  int getLayerID(TrackerHit* hit)    { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::layer]; };
  int getModuleID(TrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::module]; };
  int getSensorID(TrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::sensor]; };
  
  void setupGearGeom( const gear::GearMgr* gearMgr ) ;
  
  
  unsigned int _nLayersVTX;
  
  unsigned int _nLayersSIT;
  
  
  std::vector<float> _zLayerFTD;
  
  unsigned int _nlayersFTD;
  bool _petalBasedFTDWithOverlaps;
  int _nPhiFTD; 

  int _output_track_col_quality;
  static const int _output_track_col_quality_GOOD;
  static const int _output_track_col_quality_FAIR;
  static const int _output_track_col_quality_POOR;

  
} ;

#endif