Skip to content
Snippets Groups Projects
Commit 86039294 authored by Rosen Matev's avatar Rosen Matev :sunny:
Browse files

Merge branch 'modernize-vpoverlap' into 'master' 

Tweak VPOverlapMonitor

See merge request !2431
parents 74f695fb dfd9554b
No related branches found
No related tags found
1 merge request!2431Tweak VPOverlapMonitor
Pipeline #2654502 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp
+ 128
− 147
View file @ 86039294
/*****************************************************************************\
* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#include "DetDesc/ConditionAccessorHolder.h"
#include "Event/FitNode.h"
#include "Event/Track.h"
#include "Event/VPLightCluster.h"
#include "GaudiAlg/Consumer.h"
#include "GaudiAlg/GaudiTupleAlg.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiAlg/GaudiHistoAlg.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/Vector3DTypes.h"
#include "TrackKernel/TrackFunctors.h"
#include "VPDet/DeVP.h"
#include "fmt/format.h"
class TrackVPOverlapMonitor
: public Gaudi::Functional::Consumer<void( LHCb::Tracks const&, LHCb::VPLightClusters const&, DeVP const& ),
LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, DeVP>> {
public:
/// Standard constructor
TrackVPOverlapMonitor( const std::string& name, ISvcLocator* pSvcLocator );
void operator()( LHCb::Tracks const&, LHCb::VPLightClusters const&,
DeVP const& ) const override; ///< Algorithm execution
};
namespace {
//=============================================================================
// Calculate the residuals.
//=============================================================================
Gaudi::XYZVector getResidual( const Gaudi::XYZPoint& clusterGlobal, const LHCb::FitNode& fitnode ) {
const LHCb::State state = fitnode.unbiasedState();
return clusterGlobal - Gaudi::XYZPoint{state.x(), state.y(), state.z()};
}
// A side - left side of the VELO
bool isASide( const LHCb::FitNode& node ) {
LHCb::VPChannelID vpid = node.measurement().lhcbID().vpID();
return vpid.module() % 2 == 0;
}
} // namespace
DECLARE_COMPONENT( TrackVPOverlapMonitor )
//=============================================================================
// Standard constructor, initializes variables
//=============================================================================
TrackVPOverlapMonitor::TrackVPOverlapMonitor( const std::string& name, ISvcLocator* pSvcLocator )
: Consumer{name,
pSvcLocator,
{KeyValue{"TrackContainer", LHCb::TrackLocation::Default},
KeyValue{"ClusterContainer", LHCb::VPClusterLocation::Light},
KeyValue{"VPDetectorLocation", DeVPLocation::Default}}} {}
//=============================================================================
// Main execution
//=============================================================================
void TrackVPOverlapMonitor::operator()( const LHCb::Tracks& tracks, const LHCb::VPLightClusters& clusters,
const DeVP& det ) const {
constexpr auto nStationsInVelo = VP::NModules / 2; // 2 modules per station
for ( const LHCb::Track* track : tracks ) { // start of track-loop
const bool bwd = track->checkFlag( LHCb::Track::Flags::Backward );
const auto type = track->type();
if ( type != LHCb::Track::Types::Velo && type != LHCb::Track::Types::Long && !bwd ) { continue; }
const bool fitted = track->checkFitStatus( LHCb::Track::FitStatus::Fitted );
if ( !fitted ) continue;
// container of nodes per station (26 stations)
std::array<std::vector<const LHCb::FitNode*>, nStationsInVelo> nodescontainer{};
for ( const LHCb::FitNode* node : nodes( *track ) ) { // start cluster loop
if ( !node->hasMeasurement() ) continue;
if ( node->type() != LHCb::FitNode::Type::HitOnTrack ) continue;
node->measurement().visit(
[&]( const LHCb::Measurement::VP& vp ) { nodescontainer[vp.detElem->station()].push_back( node ); },
[]( ... ) { /* skip non-VP measurements*/ } );
namespace LHCb::Tr::Monitor {
namespace {
// Calculate the residuals.
Gaudi::XYZVector getResidual( const Gaudi::XYZPoint& clusterGlobal, const FitNode& fitnode ) {
const auto state = fitnode.unbiasedState();
return clusterGlobal - Gaudi::XYZPoint{state.x(), state.y(), state.z()};
}
// histograms per station
for ( size_t istation = 0; istation < nStationsInVelo; ++istation ) {
std::string statnum = std::to_string( istation );
if ( nodescontainer[istation].size() == 4 ) { // require 4 hits per station for the overlap track
const LHCb::FitNode* nodeA = nodescontainer[istation].front();
const LHCb::FitNode* nodeC = nodescontainer[istation].back();
if ( isASide( *nodeC ) ) { std::swap( nodeA, nodeC ); } // check if the node is A or C side
if ( isASide( *nodeA ) && !isASide( *nodeC ) ) {
const LHCb::LHCbID idA = ( nodeA->measurement() ).lhcbID();
const LHCb::LHCbID idC = ( nodeC->measurement() ).lhcbID();
auto iclusA = std::find_if( begin( clusters ), end( clusters ),
[vpx_idA = idA.vpID()]( const auto& c ) { return c.channelID() == vpx_idA; } );
if ( iclusA == end( clusters ) ) throw Error( "fail to find cluster" );
auto iclusC = std::find_if( begin( clusters ), end( clusters ),
[vpx_idC = idC.vpID()]( const auto& c ) { return c.channelID() == vpx_idC; } );
if ( iclusC == end( clusters ) ) throw Error( "fail to find cluster" );
const Gaudi::XYZPoint clusGlobalA( iclusA->x(), iclusA->y(), iclusA->z() );
const Gaudi::XYZPoint clusGlobalC( iclusC->x(), iclusC->y(), iclusC->z() );
// residuals on each module of the overlap station
double residualAX = getResidual( clusGlobalA, *nodeA ).x();
double residualAY = getResidual( clusGlobalA, *nodeA ).y();
double residualCX = getResidual( clusGlobalC, *nodeC ).x();
double residualCY = getResidual( clusGlobalC, *nodeC ).y();
// overlap residual
const double overlapResidualX = residualAX - residualCX;
const double overlapResidualY = residualAY - residualCY;
plot( residualAX, "x_residual_A_side/" + statnum + "/", "x residual A side - station " + statnum, -0.1, 0.1,
50 );
plot( residualCX, "x_residual_C_side/" + statnum + "/", "x residual C side - station " + statnum, -0.1, 0.1,
50 );
plot( residualAY, "y_residual_A_side/" + statnum + "/", "y residual A side - station " + statnum, -0.1, 0.1,
50 );
plot( residualCY, "y_residual_C_side/" + statnum + "/", "y residual C side - station " + statnum, -0.1, 0.1,
50 );
plot( overlapResidualX, "overlapresidual_X/" + statnum + "/", "x overlap residual - station " + statnum, -0.1,
0.1, 50 );
plot( overlapResidualY, "overlapresidual_Y/" + statnum + "/", "y overlap residual- station " + statnum, -0.1,
0.1, 50 );
// A (C) side - left (right) side of the VELO
enum class Side { A = 0, C = 1 };
Side side( const FitNode& node ) {
VPChannelID vpid = node.measurement().lhcbID().vpID();
return static_cast<Side>( vpid.module() % 2 );
}
} // namespace
struct VPOverlap final : Gaudi::Functional::Consumer<void( Tracks const&, VPLightClusters const&, DeVP const& ),
DetDesc::usesBaseAndConditions<GaudiHistoAlg, DeVP>> {
/// Standard constructor
VPOverlap( const std::string& name, ISvcLocator* pSvcLocator )
: Consumer{name,
pSvcLocator,
{KeyValue{"TrackContainer", TrackLocation::Default},
KeyValue{"ClusterContainer", VPClusterLocation::Light},
KeyValue{"VPDetectorLocation", DeVPLocation::Default}}} {}
///< Algorithm execution
void operator()( const Tracks& tracks, const VPLightClusters& clusters, const DeVP& det ) const override {
constexpr auto nStationsInVelo = VP::NModules / 2; // 2 modules per station
auto get_cluster_ptr = [&clusters]( const auto& node ) {
auto i =
std::find_if( begin( clusters ), end( clusters ),
[id = node.measurement().lhcbID().vpID()]( const auto& c ) { return c.channelID() == id; } );
return i != end( clusters ) ? &*i : nullptr;
};
for ( const Track* track : tracks ) { // start of track-loop
const bool bwd = track->checkFlag( Track::Flags::Backward );
const auto type = track->type();
if ( type != Track::Types::Velo && type != Track::Types::Long && !bwd ) { continue; }
const bool fitted = track->checkFitStatus( Track::FitStatus::Fitted );
if ( !fitted ) continue;
// container of nodes per station (26 stations)
std::array<std::vector<const FitNode*>, nStationsInVelo> nodescontainer{};
for ( const FitNode* node : nodes( *track ) ) { // start cluster loop
if ( !node->hasMeasurement() ) continue;
if ( node->type() != FitNode::Type::HitOnTrack ) continue;
node->measurement().visit(
[&]( const Measurement::VP& vp ) { nodescontainer[vp.detElem->station()].push_back( node ); },
[]( ... ) { /* skip non-VP measurements */ } );
}
// now plot the node and cluster positions for the hits in overlap tracks
for ( size_t inode = 0; inode < 4; ++inode ) {
plot2D( nodescontainer[istation][inode]->state().position().x(),
nodescontainer[istation][inode]->state().position().y(), 141, "x vs y node global position", -60.0,
60.0, -60.0, 60.0, 100, 100 );
plot( nodescontainer[istation][inode]->state().position().z(), 142, "z node global position", -400.0, 900.0,
50 );
const LHCb::LHCbID id = ( nodescontainer[istation][inode]->measurement() ).lhcbID();
auto cluster = std::find_if( begin( clusters ), end( clusters ),
[vp_id = id.vpID()]( const auto& c ) { return c.channelID() == vp_id; } );
plot2D( cluster->x(), cluster->y(), 143, "x vs y cluster global position ", -60.0, 60.0, -60.0, 60.0, 100,
100 );
plot( cluster->z(), 144, "z cluster global position", -400.0, 900.0, 50 );
const DeVPSensor& sensor = det.sensor( id.vpID() );
unsigned int module = sensor.module();
unsigned int station = sensor.station();
plot( module, 145, "module ", 0.0, 52.0, 52 );
plot( station, 146, "station ", 0.0, 27.0, 27 );
const Gaudi::XYZPoint clusGlobal( cluster->x(), cluster->y(), cluster->z() );
const auto clusLocal = sensor.globalToLocal( clusGlobal );
plot( clusLocal.z(), 147, "z cluster local sensor position", -10.0, 10.0, 50 );
plot2D( clusLocal.x(), clusLocal.y(), 148, "x vs y cluster local position", -10.0, 50.0, -10.0, 50.0, 100,
100 );
// histograms per station
for ( const auto& [istation, nodes] : range::enumerate( nodescontainer ) ) {
plot2D( istation, nodes.size(), 139, "#nodes vs. station", -0.5, nStationsInVelo - 0.5, -0.5, 8.5,
nStationsInVelo, 9 );
if ( nodes.size() != 4 ) continue; // require 4 hits per station for the overlap track
if ( auto nodeA = nodes.front(), nodeC = nodes.back(); side( *nodeA ) != side( *nodeC ) ) {
if ( side( *nodeC ) == Side::A ) std::swap( nodeA, nodeC );
auto iclusA = get_cluster_ptr( *nodeA );
if ( !iclusA ) throw Error( "failed to find cluster" );
auto iclusC = get_cluster_ptr( *nodeC );
if ( !iclusC ) throw Error( "failed to find cluster" );
const Gaudi::XYZPoint clusGlobalA( iclusA->x(), iclusA->y(), iclusA->z() );
const Gaudi::XYZPoint clusGlobalC( iclusC->x(), iclusC->y(), iclusC->z() );
// residuals on each module of the overlap station
double residualAX = getResidual( clusGlobalA, *nodeA ).x();
double residualAY = getResidual( clusGlobalA, *nodeA ).y();
double residualCX = getResidual( clusGlobalC, *nodeC ).x();
double residualCY = getResidual( clusGlobalC, *nodeC ).y();
// overlap residual
const double overlapResidualX = residualAX - residualCX;
const double overlapResidualY = residualAY - residualCY;
// TODO: instead of making 6 * 26 1D histograms , why not make 6 2D histograms?
plot( residualAX, fmt::format( "x_residual_A_side/{}/", istation ),
fmt::format( "x residual A side - station {}", istation ), -0.1, 0.1, 50 );
plot( residualCX, fmt::format( "x_residual_C_side/{}/", istation ),
fmt::format( "x residual C side - station {}", istation ), -0.1, 0.1, 50 );
plot( residualAY, fmt::format( "y_residual_A_side/{}/", istation ),
fmt::format( "y residual A side - station {}", istation ), -0.1, 0.1, 50 );
plot( residualCY, fmt::format( "y_residual_C_side/{}/", istation ),
fmt::format( "y residual C side - station {}", istation ), -0.1, 0.1, 50 );
plot( overlapResidualX, fmt::format( "overlapresidual_X/{}/", istation ),
fmt::format( "x overlap residual - station {}", istation ), -0.1, 0.1, 50 );
plot( overlapResidualY, fmt::format( "overlapresidual_Y/{}/", istation ),
fmt::format( "y overlap residual- station {}", istation ), -0.1, 0.1, 50 );
}
// now plot the node and cluster positions for the hits in overlap tracks
for ( const auto* inode : nodes ) {
plot2D( inode->state().position().x(), inode->state().position().y(), 141, "x vs y node global position",
-60.0, 60.0, -60.0, 60.0, 100, 100 );
plot( inode->state().position().z(), 142, "z node global position", -400.0, 900.0, 50 );
auto cluster = get_cluster_ptr( *inode );
if ( !cluster ) throw Error( "failed to find cluster" );
plot2D( cluster->x(), cluster->y(), 143, "x vs y cluster global position ", -60.0, 60.0, -60.0, 60.0, 100,
100 );
plot( cluster->z(), 144, "z cluster global position", -400.0, 900.0, 50 );
const DeVPSensor& sensor = det.sensor( inode->measurement().lhcbID().vpID() );
plot( sensor.module(), 145, "module ", 0.0, 52.0, 52 );
plot( sensor.station(), 146, "station ", 0.0, 27.0, 27 );
const auto clusLocal = sensor.globalToLocal( Gaudi::XYZPoint{cluster->x(), cluster->y(), cluster->z()} );
plot( clusLocal.z(), 147, "z cluster local sensor position", -10.0, 10.0, 50 );
plot2D( clusLocal.x(), clusLocal.y(), 148, "x vs y cluster local position", -10.0, 50.0, -10.0, 50.0, 100,
100 );
}
}
}
}
}
}
};
DECLARE_COMPONENT_WITH_ID( VPOverlap, "TrackVPOverlapMonitor" )
} // namespace LHCb::Tr::Monitor
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment