SiliconTrackingAlg.cpp~

#include "SiliconTrackingAlg.h"
#include "GearSvc/IGearSvc.h"
#include "EventSeeder/IEventSeeder.h"
#include "TrackSystemSvc/ITrackSystemSvc.h"
#include "edm4hep/MCParticle.h"
#include "edm4hep/TrackerHit.h"
//#include "edm4hep/TrackerHitPlane.h"
#include "edm4hep/Track.h"
#include "edm4hep/TrackState.h"

//#include "DataHelper/ClusterExtended.h"
//#include "DataHelper/TrackExtended.h"
//#include "DataHelper/TrackerHitExtended.h"
//#include "DataHelper/HelixClass.h"

#include <iostream>
#include <algorithm>
#include <cmath>
#include <climits>

#include <gear/GEAR.h>
#include <gear/GearMgr.h>
#include <gear/GearParameters.h>
#include <gear/VXDLayerLayout.h>
#include <gear/VXDParameters.h>
#include "gear/FTDLayerLayout.h"
#include "gear/FTDParameters.h"

#include <gear/BField.h>

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

#include "TrackSystemSvc/MarlinTrkUtils.h"
#include "TrackSystemSvc/HelixTrack.h"
#include "TrackSystemSvc/HelixFit.h"
#include "TrackSystemSvc/IMarlinTrack.h"

//#include "TrackSystemSvc/MarlinTrkDiagnostics.h"
//#ifdef MARLINTRK_DIAGNOSTICS_ON
//#include "TrackSystemSvc/DiagnosticsController.h"
//#endif

//#include "MarlinCED.h"

//#include "marlin/AIDAProcessor.h"

//---- ROOT -----
#include "TH1F.h"
#include "TH2F.h"

using namespace edm4hep ;
//using namespace marlin ;
using namespace MarlinTrk ;

using std::min;
using std::max;
using std::abs;

const int SiliconTrackingAlg::_output_track_col_quality_GOOD = 1;
const int SiliconTrackingAlg::_output_track_col_quality_FAIR = 2;
const int SiliconTrackingAlg::_output_track_col_quality_POOR = 3;

const double SiliconTrackingAlg::TWOPI = 2*M_PI;

DECLARE_COMPONENT( SiliconTrackingAlg )

SiliconTrackingAlg::SiliconTrackingAlg(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
  
  //_description = "Pattern recognition in silicon trackers";
  
  _fastfitter = new MarlinTrk::HelixFit();
  
  _encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
  
  _petalBasedFTDWithOverlaps = false;
  
  // zero triplet counters
  _ntriplets = _ntriplets_good = _ntriplets_2MCP = _ntriplets_3MCP = _ntriplets_1MCP_Bad = _ntriplets_bad = 0;

  // Input Collections
  // ^^^^^^^^^^^^^^^^^
  declareProperty("HeaderCol", _headerColHdl);
  declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
  declareProperty("VTXHitCollection", _inVTXColHdl, "Handle of the Input VTX TrackerHits collection");
  declareProperty("FTDPixelHitCollection", _inFTDPixelColHdl, "Handle of the Input FTD TrackerHits collection");
  declareProperty("FTDSpacePointCollection", _inFTDSpacePointColHdl, "Handle of the Input FTD SpacePoints collection");
  declareProperty("SITHitCollection", _inSITColHdl, "Handle of the Input SIT TrackerHits collection");
    
  // Output Collections
  // ^^^^^^^^^^^^^^^^^^
  declareProperty("SiTrackCollection", _outColHdl, "Handle of the SiTrack output collection");
  //declareProperty("TrkHitRelCollection", _outRelColHdl, "Handle of TrackerHit Track relation collection");
  // Steering parameters
  // ^^^^^^^^^^^^^^^^^^^
  
  _output_track_col_quality = _output_track_col_quality_GOOD;
  
}



StatusCode  SiliconTrackingAlg::initialize() { 
  
  _nRun = -1 ;
  _nEvt = 0 ;
  //printParameters() ;
  
  // set up the geometery needed by KalTest
  //FIXME: for now do KalTest only - make this a steering parameter to use other fitters
  auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
  if ( !_trackSystemSvc ) {
    error() << "Failed to find TrackSystemSvc ..." << endmsg;
    return StatusCode::FAILURE;
  }
  _trksystem =  _trackSystemSvc->getTrackSystem();
  
  if( _trksystem == 0 ){
    error() << "Cannot initialize MarlinTrkSystem of Type: KalTest" <<endmsg;
    return StatusCode::FAILURE;
  }
  
  _trksystem->setOption( IMarlinTrkSystem::CFG::useQMS,        _MSOn ) ;
  _trksystem->setOption( IMarlinTrkSystem::CFG::usedEdx,       _ElossOn) ;
  _trksystem->setOption( IMarlinTrkSystem::CFG::useSmoothing,  _SmoothOn) ;
  _trksystem->init() ;  
  std::cout << "fucd ==============" << _trksystem << std::endl;
#ifdef MARLINTRK_DIAGNOSTICS_ON
  
  void * dcv = _trksystem->getDiagnositicsPointer();
  DiagnosticsController* dc = static_cast<DiagnosticsController*>(dcv);
  dc->init(_MarlinTrkDiagnosticsName,_MarlinTrkDiagnosticsName, _runMarlinTrkDiagnostics);
  
#endif
  
  if(setupGearGeom()==StatusCode::FAILURE) return StatusCode::FAILURE;
  
  if (_useSIT == 0)
    _nLayers = _nLayersVTX;
  else 
    _nLayers = _nLayersVTX + _nLayersSIT;
  
  // initialise the container to have separate vectors for up to _nHitsChi2 hits.
  _tracksWithNHitsContainer.resize(_nHitsChi2);
  
  _dPhi = TWOPI/_nDivisionsInPhi;
  _dTheta = 2.0/_nDivisionsInTheta;
  _dPhiFTD = TWOPI/_nPhiFTD;
  // I leave this for the moment, but 0.3 is c/1e9.
  // For the cut it does not make too much of a difference
  double cutOnR = _cutOnPt/(0.3*_bField);
  cutOnR = 1000.*cutOnR;
  _cutOnOmega = 1/cutOnR;
  
  _output_track_col_quality = 0;
  
  return GaudiAlgorithm::initialize();
}

StatusCode SiliconTrackingAlg::execute(){ 
  
  //_current_event = evt;
  //_allHits.reserve(1000);

  _output_track_col_quality = _output_track_col_quality_GOOD;
  
  // zero triplet counters
  _ntriplets = _ntriplets_good = _ntriplets_2MCP = _ntriplets_3MCP = _ntriplets_1MCP_Bad = _ntriplets_bad = 0;
  
  // Clearing the working containers from the previous event
  // FIXME: partly done at the end of the event, in CleanUp. Make it consistent.
  //_tracksWithNHitsContainer.clear();
  //_trackImplVec.clear();
  
  //_colTrackerHits.clear();
  //_colNamesTrackerHits.clear();
  //auto header = _headerColHdl.get()->at(0);
  //int evtNo = header.getEventNumber();
  //int runNo = header.getRunNumber();
  //debug() << "Processing Run[" << runNo << "]::Event[" << evtNo << "]" << endmsg;

  _trackImplVec.reserve(100);
  _allHits.reserve(1000);

  int successVTX = InitialiseVTX();
  int successFTD = 0;
  //int successFTD = InitialiseFTD();
  if (successVTX == 1) {
    
    debug() << "      phi          theta        layer      nh o :   m :   i  :: o*m*i " << endmsg; 
    
    for (int iPhi=0; iPhi<_nDivisionsInPhi; ++iPhi) { 
      for (int iTheta=0; iTheta<_nDivisionsInTheta;++iTheta) {
        ProcessOneSector(iPhi,iTheta); // Process one VXD sector     
      }
    }
    
    debug() << "End of Processing VXD and SIT sectors" << endmsg;