TrackAnalysis.cxx

Go to the documentation of this file.

 
#include "TrigInDetAnalysisExample/ChainString.h"
 
#include "xAODEventInfo/EventInfo.h"
 
#include "TrackAnalysisStore.h"
 
#include "MsgMacros.h"
 
#include "TrackView.h"
#include "VertexView.h"
 
 
#include "TrackAnalysis.h"
#include "TBestMatch.h"
 
#include "TrackRoiSelector.h"
 
#include "TagNProbe.h"
 
#include <algorithm>
#include <limits>
#include <cmath>
#include <utility>
#include <cstdlib> 
 

IDTPM2::TrackAnalysis::TrackAnalysis( const std::string& type, const std::string& name, const IInterface* parent ) :
  AthAlgTool( type, name, parent ),
  m_tdt("Trig::TrigDecisionTool/TrigDecisionTool"),
  m_analysis(std::make_unique<AnalysisR4>(name)) 
{ }
 
 
 
 

IDTPM2::TrackAnalysis::~TrackAnalysis() = default;
 

StatusCode IDTPM2::TrackAnalysis::initialize() {
 
  MY_MSG_INFO <<  "TrackAnalysis::initialize() name: " << name() << endmsg;
 
  std::string selectChain = name().substr(name().find(".")+1,name().size());
 
  if ( selectChain.find(":")==std::string::npos ) selectChain = ":"+selectChain;
  
  ChainString chainName = selectChain; //chainName.head();

  if ( chainName.head()!="" )     m_triggerchain = chainName.head();
  if ( chainName.tail()!="" )     m_testTracks   = chainName.tail();
  if ( chainName.roi()!="" )      m_rois         = chainName.roi();
  if ( chainName.vtx()!="" )      m_vertices     = chainName.vtx();
  if ( chainName.element()!="" )  m_leg          = std::atoi(chainName.element().c_str());
  if ( chainName.extra()!="" )    m_extra        = chainName.extra();

  m_refTracks   = m_reference;
 
  m_refVertices = m_referenceVertices;

 
  MY_MSG_INFO << "analysis:   "  << name() << endmsg;
  MY_MSG_INFO << "trigger:    "  << m_triggerchain << endmsg;
  MY_MSG_INFO << "reference:  "  << m_refTracks << ( m_refTracks_tag.empty() ? "" : "\t"+ m_refTracks_tag ) << endmsg;
  MY_MSG_INFO << "test:       "  << m_testTracks << endmsg;
  if ( !m_vertices.empty() ) { MY_MSG_INFO << "vtx:       " << name() << endmsg; }
  MY_MSG_INFO << endmsg;
 
  if ( !m_tagChain.empty() ) {
    
    ChainString tagchain = m_tagChain.value();
 
    if ( m_reference_tag.empty() ) m_refTracks_tag = m_refTracks;
    else                           m_refTracks_tag = m_reference_tag;
        
    MY_MSG_INFO << "Tag and probe: " << tagchain << endmsg;
    
    MY_MSG_INFO << "Tag chain: " << tagchain.head()
        << "\t" << tagchain.roi()
        << "\t" << tagchain.element()
        << "\t" << tagchain.extra()
        << "\treferences: " << m_reference << "\t" << m_reference_tag 
        << endmsg;
 
    m_triggerchain_tag = tagchain.head();
    m_roi_tag          = tagchain.roi();
    m_extra_tag        = tagchain.extra();
    if ( !tagchain.element().empty() ) m_leg_tag = std::atoi(tagchain.element().c_str());
 
    using PtPdgHelper = PtPdgIdHelper<TrackView,xAOD::TruthParticle>;
    
    if ( m_refPdgId_tag.value() ) m_tnpPredicate = std::make_unique<PtPdgHelper>( m_minTagPT.value(), m_refPdgId_tag.value() ); 
    else                          m_tnpPredicate = std::make_unique<PtHelper<TrackView>>( m_minTagPT.value() ); 
 
    MY_MSG_INFO << "TagNProbe: " << m_triggerchain << " -> " << m_triggerchain_tag << endmsg;   
    
  }
 
  
  MY_MSG_DEBUG << "TrackAnalysis::initialize()"
           << " trigger="    << m_triggerchain
           << " testTracks=" << m_testTracks
           << " refTracks="  << m_refTracks
           << " rois="       << m_rois << endmsg;
  
  ATH_CHECK( m_eventInfoName.initialize() );
 
  ATH_CHECK( m_beamKeys.initialize() );
 
  ATH_CHECK( m_refSelector.retrieve() );
  ATH_CHECK( m_testSelector.retrieve() );
 
  //  ATH_CHECK( m_vertexQualitySelectionTool.retrieve() );
 
  ATH_CHECK( m_trackMatchingTool.retrieve() );
 
#if 0
  std::unique_ptr<TrackAnalysisDefinition> config = std::make_unique<TrackAnalysisDefinition>();
 
  config->setTestCollection( m_testTracks );
  config->setReferenceCollection( m_refTracks );
 
  config->setTestType( "Trigger" );
  config->setReferenceType( "Offline" );

  config->setUseTrigger( true );
  config->setUseOffline( true );
  config->setDoTrigNavigation( true );

  m_trkAnaDef = std::move( config );
 
#endif 
  
  //  m_trackQualitySelectionTool->setAnaDef(m_trkAnaDef.get());
 
  MY_MSG_DEBUG << "TrackAnalysis::initialize() exitting" << endmsg;
 
  m_analysis->set_monTool( &m_tool );
 
  for ( size_t i=m_fillers.size() ; i-- ; ) m_analysis->addFillTool( m_fillers[i] );
      
  m_analysis->initialise();
 
  m_initialized = true;
  
  return StatusCode::SUCCESS;
 
}
 
 

StatusCode IDTPM2::TrackAnalysis::bookHistograms()
{
  return StatusCode::SUCCESS;
}
 

StatusCode IDTPM2::TrackAnalysis::fillHistograms( const EventContext& context ) const { 
  if ( execute( context ) )  return StatusCode::SUCCESS;
  return StatusCode::FAILURE;
}
 
 
 
 
void printtracks( const std::vector<const TrackView*>& ars, const std::string& sht=""  ) {
  std::cout << "\nSUTT printtracks: " << sht << "\tsize: " << ars.size() << std::endl;
  for ( size_t i=0 ; i<ars.size() && i<10 ; i++ ) std::cout << "\t" << i << "\t" << *ars[i] << std::endl;
  std::cout << std::endl;
}
 
 
template<typename T>
std::ostream& operator<<( std::ostream& s, std::vector<T>& v ) {
  for ( size_t i=0 ; i<v.size() ; i++ ) s << "\t" << v[i];
  return s;
}
 
 
#include <fenv.h>
 
bool IDTPM2::TrackAnalysis::execute( const EventContext& context ) const {
 
  feclearexcept(FE_ALL_EXCEPT);
  fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
 
  MY_MSG_DEBUG << "Filling hists " << name() << "\ttrigger: " << m_triggerchain << " ..." << endmsg;
  
  unsigned decisionType = TrigDefs::Physics; // TrigDefs::includeFailedDecisions;
 
  if( !m_triggerchain.empty() && !m_tdt->isPassed( m_triggerchain, decisionType ) ) return true;
 
  MY_MSG_DEBUG <<  "␛[91;1m PROCESSING ------------------------------------------------------------------------------------␛[m" << endmsg;
  
  MY_MSG_DEBUG << "TrackAnalysis::execute() " << name()
           << " refTracks:     " << m_refTracks
           << "\t testTracks:  " << m_testTracks
           << "\ttriggerchain: " << m_triggerchain << endmsg;
  
 
#if 0                                                               
  if ( m_testTracks.find("HLT_")==std::string::npos && m_testTracks.find("EF_")==std::string::npos ) {
    
    const std::vector<std::string> configuredChains = m_tdt->getListOfTriggers("HLT_.*, EF_.*");
 
    for ( unsigned i=0 ; i<configuredChains.size() ; i++ ) {
      if( m_tdt->isPassed( configuredChains[i], decisionType ) ) std::cout << i << "\t" << "␛[91;1m" << "Chain " << configuredChains[i] << " pass   (ACN)␛[m" << std::endl;
      if ( configuredChains[i] == m_triggerchain ) std::cout << "Mutha ..." << configuredChains[i] << " " <<  m_triggerchain << "\n\n" << std::endl;
    }
    
    MY_MSG_INFO <<  "␛[91;1m" << configuredChains.size() << " Configured Chains" << "␛[m"
         <<"\ttrigger chain: " << m_triggerchain << " :: "
         << m_tdt->isPassed( m_triggerchain, decisionType ) << endmsg;
  }
#endif
  
  IDTPM2::TrackAnalysisStore trackAnaStore( "_duff", context );
 
  if ( !m_useDecorations ) trackAnaStore.setBeamKeys( &m_beamKeys );
    
  //  trackAnaStore.setTDT( &m_tdt );
 
  RoiSelector roiSelector(&m_tdt);
 
  trackAnaStore.fillViewRef( m_refTracks );
 
  //  printtracks( trackAnaStore.reference_tracks(), "\ttrackAnaStore: reference" );  
  
 
  bool isTrigger = m_testTracks.find("HLT_")!=std::string::npos || m_testTracks.find("EF_")!=std::string::npos; 
  
  if ( !isTrigger ) trackAnaStore.fillViewTest( m_testTracks );

  std::vector<TrigCompositeUtils::LinkInfo<TrigRoiDescriptorCollection> > selectedRois; 

  std::vector<TrigCompositeUtils::LinkInfo<TrigRoiDescriptorCollection> > tagRois;  
  
  size_t selectedRoisSize = 1; // by default only one "dummy" RoI, i.e. for offline analysis 

  
  TrackRoiSelector roiFilter;
 
  std::unique_ptr<TagNProbe> tnp = 0;

 
  if ( !m_tagChain.empty() ) tnp = std::make_unique<TagNProbe>( m_minTagMass, m_maxTagMass );

  if ( !m_triggerchain.empty() ) { 
    selectedRois     = roiSelector.retrieveRois( m_triggerchain, m_rois ); 
    selectedRoisSize = selectedRois.size();
  }


  std::unique_ptr<TrackAnalysisStore> tagStore = 0;
 
  if ( tnp ) {


      const std::vector<const TrackView*>*  vref_ptr;
 
      MY_MSG_DEBUG << "reftracks: " << m_refTracks << "\t:" << m_refTracks_tag << endmsg;
      
      if ( m_refTracks!=m_refTracks_tag ) {
        tagStore = std::make_unique<TrackAnalysisStore>( "tagstore", context );
    tagStore->fillViewRef( m_refTracks_tag );
    vref_ptr = &tagStore->reference_tracks();
      }
      else {
        vref_ptr = &trackAnaStore.reference_tracks();
      }
 
      const std::vector<const TrackView*>& vref = *vref_ptr;
      
      tagRois = roiSelector.retrieveRois( m_triggerchain_tag, m_roi_tag, m_leg_tag );
 
      for ( auto roi : tagRois ) {
        ElementLink<TrigRoiDescriptorCollection> roi_link = roi.link;
        if( ! roi_link.isValid() ) continue;
    std::vector<const TrackView*> v = roiFilter.getTracks( vref, *roi_link.cptr() );
    tnp->append_tags( v.begin(), v.end(), *m_tnpPredicate.get() );
      }
 
      MY_MSG_DEBUG << "tnp size() : " << tnp->size() << "\t" << vref.size() << endmsg;
 
#if 0
      const std::vector<const TrackView*>& v = trackAnaStore.reference_tracks();
      for ( size_t i=0 ; i<v.size() ; i++ ) {
    bool ok = tnp->select( v[i] );
    std::vector<double> m = tnp->masses( v[i] );
    std::cout << "\ttnp: " << ok << "\tm: " << m.size() << m << "\t " << *v[i] << std::endl;
      }
#endif  
 
  }
  
 
  TrackAnalysisStore      roiStore(      "roi", &m_tdt, context );  
  TrackAnalysisStore selectedStore( "selected", &m_tdt, context );

  if ( !m_vertices.empty() && !m_refVertices.empty() ) selectedStore.fillVertexViewRef( m_refVertices );

 
  //  std::cout << "TrackAnalysis: ref vtx size: " << selectedStore.refvertices().size() << std::endl;
 
  //  const std::vector<const VertexView*> vtx = selectedStore.refvertices();
  
  //  for ( int iv=vtx.size() ; iv-- ; ) {
  //   if ( vtx[iv]->vertexType()!=0 ) std::cout << "\t" << iv << "\tvtx[" << iv << "]: " << *vtx[iv] << std::endl;
  //  }
      
 
  
  if ( selectedStore.fillEventInfo( m_eventInfoName.key() ).isFailure() ) return false;
 
  MY_MSG_DEBUG << "Processing event = " << selectedStore.eventInfo()->eventNumber() << endmsg;
 
  for( size_t ir=selectedRoisSize ; ir-- ;  ) {
 
    MY_MSG_DEBUG << "\tRoi: " << ir << endmsg; 
    
    if( !m_triggerchain.empty() ) { 

        selectedStore.roiIndex(ir);
      
        ElementLink<TrigRoiDescriptorCollection> roi_link = selectedRois[ir].link;

        if( ! roi_link.isValid() ) {
          MY_MSG_WARNING << "Found non-valid RoI ElementLink" << endmsg;
          continue;
        }
 
    const TrigRoiDescriptor* const* roiptr = roi_link.cptr();
 
    if ( !roiptr ) {
      MY_MSG_WARNING << "could not fetch RoiDescriptor: " << m_rois << endmsg;
      continue;
    }

    if ( isTrigger ) { 
 
      roiStore.fillViewTest( m_testTracks, roi_link  );

      if ( !m_vertices.empty() ) selectedStore.fillVertexViewTest( m_vertices, roi_link );

      //      const std::vector<const VertexView*> vtx = selectedStore.testvertices();
  
      //      for ( int iv=vtx.size() ; iv-- ; ) {
      //         std::cout << "\t\t" << iv << "\tvtx[" << iv << "]: " << *vtx[iv] << std::endl;
      //      }
      
      //      printtracks( roiStore.test_tracks(),            "\nroiStore: test tracks" );
      //      if ( !m_vertices.empty() ) fillroivertices( roi_link, m_vertices );
            
      //        selectedStore.fill<VertexView,xAOD::VertexContainer>( selectedStore.testvertices(), 
      //                                  roi_link, m_vertices );
 
    }

    //  printtracks( trackAnaStore.test_tracks(), "test (should be ok - roi retrieve II)" );  
      
    
        roiStore.reference_tracks() = roiFilter.getTracks( trackAnaStore.reference_tracks(), *roiptr );

 
    }
    else {
      roiStore.reference_tracks() = trackAnaStore.reference_tracks();
      roiStore.test_tracks()      = trackAnaStore.test_tracks();
    }
 
    //    printtracks( roiStore.reference_tracks(), "\troiStore: reference" );  
    //    printtracks( roiStore.test_tracks(), "\troiStore: test" );
    
    MY_MSG_DEBUG << "SUTT: TrackSelection ..." << endmsg;
 
    selectedStore.reference_tracks() = m_refSelector->selectTracks( roiStore.reference_tracks(), tnp.get() );
 
    //    printtracks( selectedStore.reference_tracks(), "\tselectedStore: reference" );  
    
    selectedStore.test_tracks() = m_testSelector->selectTracks( roiStore.test_tracks() );
 
    //    printtracks( selectedStore.test_tracks(), "\tselectedStore: test" );
    
    // TDeltaRMatcher<const TrackView> matcher( "DeltaRMatcher", 0.1 ); /// FIXME: turn this hard coded parameter in to a job Property
 
    std::unique_ptr<TIDA::Associator<const TrackView> > matcher = m_trackMatchingTool->uclone();
    
    //    MY_MSG_DEBUG << "n matched: " << matcher->size() << endmsg; 

    //  matchertype mt = matcher.uclone(); 
    matcher->match( selectedStore.reference_tracks(), selectedStore.test_tracks() );
    
    // matcher.match( tref, ttest );
    
    //    MY_MSG_DEBUG <<  "SUTT: n matched: " << matcher->size() << endmsg;
    
    selectedStore.set_matcher( matcher.get() );
 
    TDeltaZMatcher<const VertexView> vertex_matcher( "DeltaZMatcher", 1.0 ); 
 
    if ( !m_vertices.empty() ) {
      vertex_matcher.match( selectedStore.refvertices(), selectedStore.testvertices() );
      selectedStore.set_vertex_matcher( &vertex_matcher );
    }
 
    if ( tnp ) selectedStore.addTagNProbe( tnp.get() );
    
    //    selectedStore.refMap( trackAnaStore );
  
 
    m_analysis->execute( selectedStore );
 
  }
 
  return true;
  
}
 

StatusCode IDTPM2::TrackAnalysis::proc() {
 
  MY_MSG_DEBUG << "Finalising plots" << endmsg;

  // if( endOfRunFlag() ) m_analysis->finalize()
 
  m_analysis->finalize();
  
  MY_MSG_DEBUG << "Successfully finalized hists" << endmsg;
 
  return StatusCode::SUCCESS;
}