CTPPSDiamondTrackRecognition.cc

Go to the documentation of this file.

/****************************************************************************
 *
 * This is a part of CTPPS offline software.
 * Authors:
 *   Laurent Forthomme (laurent.forthomme@cern.ch)
 *   Nicola Minafra (nicola.minafra@cern.ch)
 *
 ****************************************************************************/

#include "RecoCTPPS/TotemRPLocal/interface/CTPPSDiamondTrackRecognition.h"

#include <cmath>
#include <cstdio>
#include <algorithm>

//----------------------------------------------------------------------------------------------------

const std::string CTPPSDiamondTrackRecognition::pixelEfficiencyDefaultFunction_ = "(x>[0]-0.5*[1])*(x<[0]+0.5*[1])";

CTPPSDiamondTrackRecognition::CTPPSDiamondTrackRecognition( const edm::ParameterSet& iConfig ) :
  threshold_           ( iConfig.getParameter<double>( "threshold" ) ),
  thresholdFromMaximum_( iConfig.getParameter<double>( "thresholdFromMaximum" ) ),
  resolution_          ( iConfig.getParameter<double>( "resolution" ) ),
  sigma_               ( iConfig.getParameter<double>( "sigma" ) ),
  startFromX_          ( iConfig.getParameter<double>( "startFromX" ) ),
  stopAtX_             ( iConfig.getParameter<double>( "stopAtX" ) ),
  yPositionInitial_    ( iConfig.getParameter<double>( "yPosition" ) ),
  yWidthInitial_       ( iConfig.getParameter<double>( "yWidth" ) ),
  hit_f_( "hit_TF1_CTPPS", iConfig.getParameter<std::string>( "pixelEfficiencyFunction" ).c_str(), startFromX_, stopAtX_ )
{
  if ( sigma_ == 0. ) {
    hit_f_ = TF1( "hit_TF1_CTPPS", pixelEfficiencyDefaultFunction_.c_str(), startFromX_, stopAtX_ ); // simple step function
  }
}

//----------------------------------------------------------------------------------------------------

CTPPSDiamondTrackRecognition::~CTPPSDiamondTrackRecognition()
{}

//----------------------------------------------------------------------------------------------------

void
CTPPSDiamondTrackRecognition::clear()
{
  hitParametersVectorMap_.clear();
  mhMap_.clear();
  yPosition_ = yPositionInitial_;
  yWidth_ = yWidthInitial_;
  zPosition_ = zPositionInitial_;
  zWidth_ = zWidthInitial_;
}

//----------------------------------------------------------------------------------------------------

void
CTPPSDiamondTrackRecognition::addHit( const CTPPSDiamondRecHit& recHit )
{
  // store hit parameters
  hitParametersVectorMap_[recHit.getOOTIndex()].emplace_back( recHit.getX(), recHit.getXWidth() );

  // check y
  if ( yPosition_ == yPositionInitial_ && yWidth_ == yWidthInitial_ ) {
    yPosition_ = recHit.getY();
    yWidth_ = recHit.getYWidth();
  }

  // check z
  if ( zPosition_ == zPositionInitial_ && zWidth_ == zWidthInitial_ ) {
    zPosition_ = recHit.getZ();
    zWidth_ = recHit.getZWidth();
  }
}

//----------------------------------------------------------------------------------------------------

int
CTPPSDiamondTrackRecognition::produceTracks( edm::DetSet<CTPPSDiamondLocalTrack>& tracks )
{
  int number_of_tracks = 0;
  const double inv_resolution = 1./resolution_;

  for ( auto const& oot : hitParametersVectorMap_ ) {
    std::vector<float> hit_profile( ( stopAtX_-startFromX_ ) * inv_resolution, 0. );
    for ( auto const& param : oot.second ) {
      hit_f_.SetParameters( param.center, param.width, sigma_ );
      for ( unsigned int i = 0; i < hit_profile.size(); ++i ) {
        hit_profile[i] += hit_f_.Eval( startFromX_ + i*resolution_ );
      }
    }

    float maximum = 0.;
    bool below = true; // start below the threshold
    int track_start_n = 0;

    for ( unsigned int i = 0; i < hit_profile.size(); ++i ) {
      if ( below && hit_profile[i] >= threshold_ ) { // going above the threshold
        track_start_n = i;
        maximum = 0;
        below = false;
      }
      if ( !below ) {
        if ( hit_profile[i] > maximum ) {
          maximum = hit_profile[i];
        }
        if ( hit_profile[i] < threshold_ ) { // going back below the threshold
          below = true;

          // go back and use new threshold
          const float threshold = maximum - thresholdFromMaximum_;
          for ( unsigned int j = track_start_n; j <= i; ++j ) {
            if ( below && hit_profile[j] >= threshold ) { // going above the threshold
              track_start_n = j;
              below = false;
            }
            if ( !below && hit_profile[j] < threshold ) { // going back below the threshold
              below = true;

              //store track
              math::XYZPoint pos0_sigma( ( j-track_start_n )*resolution_*0.5, yWidth_ * 0.5, zWidth_ * 0.5 );
              math::XYZPoint pos0( startFromX_ + track_start_n*resolution_ + pos0_sigma.X(), yPosition_, zPosition_ );
              int mult_hits = 0;
              if ( mhMap_.find( oot.first ) != mhMap_.end() ) mult_hits = mhMap_[oot.first];

              CTPPSDiamondLocalTrack track( pos0, pos0_sigma, 0., 0., oot.first, mult_hits );
              track.setValid( true );
              tracks.push_back( track );
              ++number_of_tracks;
            }
          }
        }
      }
    }
  }

  return number_of_tracks;
}