Resolutions.cxx

Go to the documentation of this file.

 
 

 
#include "AsgMessaging/StatusCode.h"
#include "AthenaBaseComps/AthMessaging.h"
 
#include "GaudiKernel/SystemOfUnits.h"
 
#include "InDetPhysValMonitoring/ResolutionHelper.h"

#include "Resolutions.h"
#include "TrackView.h"
 
#include "MsgMacros.h"
 
#include <set>
 
 

IDTPM2::Resolutions::Resolutions( const std::string& name, ToolHandle<GenericMonitoringTool>* m ) 
  : IPlotBase( name, m ), m_method(0) 
{ }
 
 

void IDTPM2::Resolutions::initialize()
{
  MY_MSG_DEBUG <<  "Booking resolutions: " << name()  << endmsg; 
  
  for( unsigned int i=1; i<NPARAMS; i++ ) {
    
    m_pull[i]         = TIDA::Histogram<float>(   monTool(), "pull_"+m_paramName[i] );
    
    m_pull[i]         = TIDA::Histogram<float>(   monTool(), "pull_"+m_paramName[i] );
    m_res[i]          = TIDA::Histogram<float>(   monTool(), "res_"+m_paramName[i] );
    m_sigma[i]        = TIDA::Histogram<float>(   monTool(), "sigma_"+m_paramName[i] );
    m_significance[i] = TIDA::Histogram<float>(   monTool(), "significance_"+m_paramName[i] );
    m_corr[i]         = TIDA::Histogram2D<float>( monTool(), "corr_test_vs_ref_"+m_paramName[i] ); 
 
    m_resolution_vs_pt[i]  = TIDA::Histogram2D<float>( monTool(), "2D", "resolution_"+m_paramName[i]+"_vs_pt" );
    m_resolution_vs_eta[i] = TIDA::Histogram2D<float>( monTool(), "2D", "resolution_"+m_paramName[i]+"_vs_eta" );
    m_resolution_vs_z0[i]  = TIDA::Histogram2D<float>( monTool(), "@D", "resolution_"+m_paramName[i]+"_vs_z0" );
    
    for( unsigned int j=0; j<NPARAMSOUT; j++ ) {

      m_resHelper[i][j] = TIDA::Histogram<float>( monTool(), "resHelper_"+m_paramName[i]+"_vs_"+m_paramName[j] );
      m_reswidth[i][j]  = TIDA::Histogram<float>( monTool(), "resolution_"+m_paramName[i]+"_vs_"+m_paramName[j] );
      m_resmean[i][j]   = TIDA::Histogram<float>( monTool(), "resmean_"+m_paramName[i]+"_vs_"+m_paramName[j] );
      
      m_pullHelper[i][j] = TIDA::Histogram<float>( monTool(), "pullHelper_"+m_paramName[i]+"_vs_"+m_paramName[j] ); 
      m_pullmean[i][j]   = TIDA::Histogram<float>( monTool(), "pullmean_"+m_paramName[i]+"_vs_"+m_paramName[j] );
      m_pullwidth[i][j]  = TIDA::Histogram<float>( monTool(), "pullwidth_"+m_paramName[i]+"_vs_"+m_paramName[j] );
      
    } // close NPARAMSOUT loop
  } // close NPARAMS loop
 
}
 


 
namespace { 
  double deltaPhi( double dphi ) {
    while ( dphi>M_PI  ) dphi -= 2*M_PI;
    while ( dphi<-M_PI ) dphi += 2*M_PI;
    return dphi;
  }
}
 
 
void IDTPM2::Resolutions::execute( const TrackAnalysisStore& collections ) const {
  
  
  
  
  //  const std::vector<const TrackView*>& refTracks  = collections.reference_tracks();
  //  const std::vector<const TrackView*>& testTracks = collections.test_tracks();
  
  const TIDA::Associator<const TrackView>* matcher = collections.matcher();
  
  std::array<float,NPARAMS>       refP = {};
  std::array<float,NPARAMS>  refErrorP = {};
  
  std::array<float,NPARAMS>      testP = {};
  std::array<float,NPARAMS> testErrorP = {};
  
  
  
  const TIDA::Associator<const TrackView>::map_type&  matches = matcher->matched();
  
  
  TIDA::Associator<const TrackView>::map_type::const_iterator itr = matches.begin();
 
  for ( ; itr!=matches.end() ; itr++ ) {  //cppcheck-suppress postfixOperator
    
    const TrackView* ref  = itr->first;
    const TrackView* test = itr->second;
    
    //  for ( size_t i=0 ; i<refTracks.size() ; i++ ) { 
    
    
    getTrackParameters(  *ref,  refP,  refErrorP );
    getTrackParameters( *test, testP, testErrorP );

    
    for( unsigned i=1 ; i<NPARAMS ; i++ ) {
      
      float residual = testP[i] - refP[i];
 
      if( i==PHI ) residual = deltaPhi( testP[i] - refP[i] ); // angular difference
      
 
      if ( testErrorP[i]!=0 )  { 
    
    float pull = residual/testErrorP[i];
    
    float sigma_test        = testErrorP[i];
 
    float significance_test = testP[i]/testErrorP[i];
    
    m_pull[i].Fill( pull ); // , weight );

    m_significance[i].Fill( significance_test ); // , weight );
 
    if ( i==QOVERPT ) {
      residual    = testP[i]/refP[i] - 1;    // Relative q/pt resolution
      sigma_test *= 1/std::fabs( testP[i] ); // relative q/pt error
    }
 
    m_sigma[i].Fill( sigma_test ); //, weight );
 
      }
        
      //    m_corr[i].Fill( refP[i], testP[i], weight ); /// !!! 2D histograms not implemented yet - next thing to do
 
      m_res[i].Fill( residual ); //, weight );
 
      m_corr[i].Fill( refP[i], testP[i] );
 
      m_resolution_vs_pt[i].Fill(  ref->pt()*0.001, residual );
      m_resolution_vs_eta[i].Fill( ref->eta(),      residual );
      m_resolution_vs_z0[i].Fill(  ref->z0(),       residual );
      
      // for( unsigned int j=0; j<NPARAMSOUT; j++ ) {
      //      m_resHelper[i][j],Fill(  refP[j], residual ); //, weight ); 
      //      m_pullHelper[i][j].Fill( refP[j], pull );     //, weight ); 
      // }
    }
  }
}
 
 

 
void IDTPM2::Resolutions::finalize() {
 
  MY_MSG_DEBUG <<  "Finalising resolution plots: " << name() << endmsg;
 
  std::string realpath = realPath( "/Resolutions/" ); 
  
  std::string names[3] = { "pt", "eta", "z0" };
    
  for( unsigned i=1 ; i<NPARAMS ; i++ ) {
    for ( int j=0 ; j<3 ; j++ ) {
      processResolution( realpath+"resolution_"+m_paramName[i]+"_vs_"+names[j]+"/" );
    }
  }
  
#if 0
  IDPVM::ResolutionHelper resolutionHelper;
  IDPVM::ResolutionHelper::methods method = IDPVM::ResolutionHelper::methods( m_method );

  for( unsigned int i=1; i<NPARAMS; i++ ) {
    for( unsigned int j=0; j<NPARAMSOUT; j++ ) {
      
      std::cout << "\t"  <<  m_resHelper[i][j].name() << std::endl;
      std::cout << "\t"  <<  m_reswidth[i][j].name()  << std::endl;
      std::cout << "\t"  <<  m_resmean[i][j].name()   << std::endl;
      
      std::cout << "\t"  <<  m_pullHelper[i][j].name() << std::endl;
      std::cout << "\t"  <<  m_pullwidth[i][j].name()  << std::endl;
      std::cout << "\t"  <<  m_pullmean[i][j].name()   << std::endl;    
      
      //      resolutionHelper.makeResolutions(  m_resHelper[i][j], m_reswidth[i][j], m_resmean[i][j],  method );
      //      resolutionHelper.makeResolutions( m_pullHelper[i][j], m_pullwidth[i][j, m_pullmean[i][j], method );
    }
  }
#endif
  
}
 

void IDTPM2::Resolutions::getTrackParameters( const TrackView& p,
                          std::array<float,NPARAMS>& params,
                          std::array<float,NPARAMS>& errors ) const
{
 
  params.fill(0);
  
  params[ D0 ]      = p.d0();
  params[ Z0 ]      = p.z0();
  params[ QOVERP ]  = p.qOverP() * Gaudi::Units::GeV;
  params[ QOVERPT ] = p.qOverP()*std::sin(p.theta()) * Gaudi::Units::GeV;
  params[ THETA ]   = p.theta();
  params[ PHI ]     = p.phi();
  params[ PT ]      = p.pt() / Gaudi::Units::GeV;
  params[ Z0SIN ]     = p.z0()*std::sin(p.theta()); 
  params[ ETA ]       = p.eta();
  params[ INCLUSIVE ] = 1;
 
  errors.fill(0);
  
  if ( p.is<xAOD::TrackParticle>() ) { 
    errors[ D0 ]        =  p.dd0();      //error( p, Trk::d0 );
    errors[ Z0 ]        =  p.dz0();      //, Trk::z0 );
    errors[ QOVERP ]    =  p.dqOverP();   // p, Trk::qOverP ) * Gaudi::Units::GeV;
    
    // errors[ QOVERPT ]   = qOverPTError( p ) * Gaudi::Units::GeV;
    errors[ QOVERPT ]   =  0;   // for the time being
    
    errors[ THETA ]     =  p.dtheta();   // error( p, Trk::theta );
    errors[ PHI ]       =  p.dphi();     // error( p, Trk::phi );
    errors[ PT ]        =  p.dpt();      // pTError( p ) / Gaudi::Units::GeV;
    //   errors[ Z0SIN ]     = z0SinThetaError( p ); ///
    errors[ Z0SIN ]     =  0; 
    errors[ ETA ]       =  p.deta();     // etaError( p );
    errors[ INCLUSIVE ] =  0;
  }
}