FullLDCTrackingAlg.cpp

#include "FullLDCTrackingAlg.h"

#include "DataHelper/Navigation.h"
#include "Tracking/TrackingHelper.h"

#include <GearSvc/IGearSvc.h>

#include <edm4hep/TrackerHitConst.h>
#include <edm4hep/TrackerHit.h>
#include <edm4hep/TrackConst.h>
#include <edm4hep/Track.h>

#include <iostream>
#include <algorithm>
#include <memory>

#include <math.h>
#include <map>

//#include "DataHelper/ClusterShapes.h"
#include <gear/GEAR.h>
#include <gear/GearParameters.h>
#include <gear/BField.h>
#include <gear/VXDLayerLayout.h>
#include <gear/VXDParameters.h>
#include "gear/FTDLayerLayout.h"
#include "gear/FTDParameters.h"
#include <gear/TPCParameters.h>
#include <gear/PadRowLayout2D.h>

#include "TrackSystemSvc/IMarlinTrack.h"
#include "TrackSystemSvc/IMarlinTrkSystem.h"
#include "TrackSystemSvc/MarlinTrkUtils.h"
#include "TrackSystemSvc/LCIOTrackPropagators.h"

#include <UTIL/LCTOOLS.h>
//#include <UTIL/LCRelationNavigator.h>

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

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

#include <climits>
#include <cmath>

#include "gsl/gsl_randist.h"
#include "gsl/gsl_cdf.h"

#include <vector>
#include <bitset>

typedef std::vector<edm4hep::ConstTrackerHit> ConstTrackerHitVec;

using namespace edm4hep ;
using namespace MarlinTrk ;

/** debug printout helper method */
std::string toString( int iTrk, edm4hep::ConstTrack tpcTrack, float bField=3.5 ) {
  
  int   nHits    = int( tpcTrack.trackerHits_size() );
  float d0TPC    = getD0(tpcTrack);
  float z0TPC    = getZ0(tpcTrack);
  float omegaTPC = getOmega(tpcTrack);
  float phi0TPC  = getPhi(tpcTrack);
  float tanLTPC  = getTanLambda(tpcTrack);
  float Chi2TPC  = tpcTrack.getChi2()/float(tpcTrack.getNdf());
  int   ndfTPC   = tpcTrack.getNdf();

  int nlinkedTracks = tpcTrack.tracks_size();
  
  
  HelixClass helixTPC;
  
  helixTPC.Initialize_Canonical(phi0TPC,d0TPC,z0TPC,omegaTPC,tanLTPC, bField);
  
  char strg[200];
  
  float pxTPC = helixTPC.getMomentum()[0];
  float pyTPC = helixTPC.getMomentum()[1];
  float pzTPC = helixTPC.getMomentum()[2];
  const float ptot = sqrt(pxTPC*pxTPC+pyTPC*pyTPC+pzTPC*pzTPC);

  sprintf(strg,"%3i   %5i %9.3f  %9.3f  %9.3f  %7.2f  %7.2f  %7.2f %4i %4i %8.3f %8i",iTrk,tpcTrack.id(),
	  ptot, d0TPC,z0TPC,pxTPC,pyTPC,pzTPC,nHits,ndfTPC,Chi2TPC,nlinkedTracks);

  return std::string( strg ) ;
}

DECLARE_COMPONENT( FullLDCTrackingAlg )
FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {  
  // _description = "Performs full tracking in ILD detector" ;  
  
  _encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
  
  // Input tracker hit collections
  declareProperty("FTDPixelTrackerHits", _FTDPixelTrackerHitColHdl, "handler of FTD Pixel Hit Collection");
  declareProperty("FTDSpacePoints", _FTDSpacePointColHdl, "FTD FTDSpacePoint Collection");
  declareProperty("VTXTrackerHits", _VTXTrackerHitColHdl, "VTX Hit Collection");
  declareProperty("SITTrackerHits", _SITTrackerHitColHdl, "SIT Hit Collection");
  declareProperty("SETTrackerHits", _SETTrackerHitColHdl, "SET Hit Collection");
  //declareProperty("ETDTrackerHits", _ETDTrackerHitColHdl, "ETD Hit Collection");
  declareProperty("TPCTrackerHits", _TPCTrackerHitColHdl, "TPC Hit Collection");
  declareProperty("SITRawHits", _inSITRawColHdl, "SIT Raw Hit Collection of SpacePoints");
  declareProperty("SETRawHits", _inSETRawColHdl, "SET Raw Hit Collection of SpacePoints");
  declareProperty("FTDRawHits", _inFTDRawColHdl, "FTD Raw Hit Collection of SpacePoints");
  //declareProperty("VTXRawHits", _inVXDRawColHdl, "VXD SimTrackerHit collection");
  
  // Input track collections
  declareProperty("TPCTracks", _TPCTrackColHdl, "TPC Track Collection"); 
  declareProperty("SiTracks", _SiTrackColHdl, "Si Track Collection");
  
  // Input relation collections
  
  /*
  registerInputCollection(LCIO::LCRELATION,
                          "TPCTracksMCPRelColl",
                          "TPC Track to MCP Relation Collection Name",
                          _TPCTrackMCPCollName,
                          std::string("TPCTracksMCP"));
  
  registerInputCollection(LCIO::LCRELATION,
                          "SiTracksMCPRelColl",
                          "Si Track to Collection",
                          _SiTrackMCPCollName,
                          std::string("SiTracksMCP"));
  */ 
  // Output track collection
  declareProperty("OutputTracks", _OutputTrackColHdl, "Full LDC track collection name");
}



StatusCode FullLDCTrackingAlg::initialize() { 
  
  // usually a good idea to
  // printParameters();  
  _nRun = -1 ;
  _nEvt = 0 ;
  PI = acos(-1.);
  PIOVER2 = 0.5*PI;
  TWOPI = 2*PI;
  
  // 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() << "Fail to find TrackSystemSvc ..." << endmsg;
  } 
  
  _trksystem = _trackSystemSvc->getTrackSystem(this);
  if( _trksystem == 0 ){
    error() << "Cannot initialize MarlinTrkSystem of Type: KalTest" <<endmsg;
    return StatusCode::FAILURE;
  }

  _trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS,        _MSOn ) ;
  _trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::usedEdx,       _ElossOn) ;
  _trksystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing,  _SmoothOn) ;
  _trksystem->init() ;
    
  this->setupGearGeom();
  return GaudiAlgorithm::initialize();
}

/*
void FullLDCTrackingAlg::processRunHeader( LCRunHeader* run) { 
  
  _nRun++ ;
  _nEvt = 0;
  streamlog_out(DEBUG5) << endmsg;
  streamlog_out(DEBUG5) << "FullLDCTrackingAlg ---> new run : run number = " << run->getRunNumber() << endmsg;
  
} 
*/

StatusCode FullLDCTrackingAlg::execute() { 
  
  // debug() << endmsg;
  debug() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << "  event = " << _nEvt << endmsg;
  // debug() << endmsg;
  auto outCol = _OutputTrackColHdl.createAndPut();
  
  prepareVectors();
  debug() << "************************************PrepareVectors done..." << endmsg;

  debug() << "************************************Merge TPC/Si ..." << endmsg;
  
  MergeTPCandSiTracks();
  debug() << "************************************Merging done ..." << endmsg;

  MergeTPCandSiTracksII();
  debug() << "************************************Merging II done ..." << endmsg;

  Sorting(_allCombinedTracks);
  debug() << "************************************Sorting by Chi2/NDF done ..." << endmsg;