Commit 24db9f99 authored by Niels Van Eldik's avatar Niels Van Eldik Committed by Graeme Stewart
Browse files

small bug fix (MuonSegmentCleaner-01-00-00)

parent ff52fbc4
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef MuonSegmentCleaner_MuonPhiHitSelector_H
#define MuonSegmentCleaner_MuonPhiHitSelector_H
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgTool.h"
#include "MuonRecToolInterfaces/IMuonHitSelector.h"
#include "MuonRecToolInterfaces/IMuonClusterOnTrackCreator.h"
#include "MuonRecToolInterfaces/IMuonCompetingClustersOnTrackCreator.h"
#include "GaudiKernel/ToolHandle.h"
class StoreGate;
class RpcIdHelper;
class MdtIdHelper;
class CscIdHelper;
class TgcIdHelper;
class StoreGateSvc;
class Identifier;
namespace MuonGM {
class MuonDetectorManager;
}
namespace Trk {
class PrepRawData;
class RIO_OnTrack;
}
class MuonPhiHitSelector : public AlgTool, virtual public Muon::IMuonHitSelector
{
public:
/** constructor */
MuonPhiHitSelector(const std::string&,const std::string&,const IInterface*);
/** destructor */
virtual ~MuonPhiHitSelector();
/** to initiate private members */
virtual StatusCode initialize();
/** to delete private members */
virtual StatusCode finalize();
/** @brief Selects and builds a cleaned vector of RIO
fits the associatedHits and build new RIOs, if m_competingRios true then for ambiguous hits competing rios are built
*/
virtual std::vector<const Trk::MeasurementBase*>* select_rio( const double pmom, const std::vector<const Trk::RIO_OnTrack*>& associatedHits,
const std::vector<const Trk::PrepRawData*>& unassociatedHits ) const;
/** return fitted phi */
virtual double getPhi()const;
private:
StoreGateSvc* m_storeGate; //!< Pointer to storegate
const MuonGM::MuonDetectorManager* m_detMgr; //!< Pointer to detector manager
const RpcIdHelper* m_rpcIdHelper; //!< Pointer to RPC id helper
const TgcIdHelper* m_tgcIdHelper; //!< Pointer to TGC id helper
const CscIdHelper* m_cscIdHelper; //!< Pointer to CSC id helper
const MdtIdHelper* m_mdtIdHelper; //!< Pointer to MDT id helper
/** Toolhandle to CompetingRIOsOnTrackTool creator */
ToolHandle<Muon::IMuonCompetingClustersOnTrackCreator> m_competingRIOsOnTrackTool;
/** Toolhandle to ClusterOnTrackTool creator */
ToolHandle<Muon::IMuonClusterOnTrackCreator> m_clusterCreator;
/** flag to print out debugging information */
bool m_debug;
/** flag to print out a summary of what comes in and what comes out */
bool m_summary;
/** flag for use of cosmics, straight line model will be used, no interaction point constraint */
bool m_cosmics;
/** flag that performs a clusterization and return clusters (default: false) */
bool m_makeClusters;
/** flag that build competing rios on track for amibguous trigger hits (default: false) */
bool m_competingRios;
/** fitted phi value */
mutable double m_phi;
/** fit method curved track model */
void fitRecPhi( const double pmom, const std::vector<Identifier> & phiId, const std::vector<double> & phiHitx, const std::vector<double> & phiHity, const std::vector<double> & phiHitz, const std::vector<double> & phiError, std::vector<int> & quality, const int nphi, std::vector<double> & phiPull, std::vector<int> & phiMult, std::vector<int> & phiSelect, double & chi2, double & r0, double & phi, std::vector<double> & errorM, int & nfit) const;
/** fit method straight line model */
void fitPhiSL(const double pmom, const std::vector<Identifier> & id, const std::vector<double> & hitx, const std::vector<double> & hity, const std::vector<double> & hitz, const std::vector<double> & error, std::vector<int> & select, const int n, std::vector<double> & pull, int & imax , double & chi2, double & r0, double & phi , std::vector<double> & errorM, bool fast) const;
/** @brief clusterization method
Use hits (select > 0) and pulls to make clusters
Inputs
id = identifiers hits
hitx hity hitz = position in space
error = associated error (in x-y plane)
pull (from fit)= residual (hit -fit) /error
select = > 0 for selected hits
n = total number of hits
fast = true = fast fit without scattering centres and no error propagation
false = fit with scattering centres and error propagation
Outputs
clusterX Y Z = cluster positions
clusterError = errors
clusterId = cluster identifier (smallest pull)
clusterHits = number of hits per cluster
ncl = number of clusters
chi2 = total chi2
r0 = perigee parameter of fit (0,0)
phi = azimuthal angle of fit at perigee
*/
void clusterPhi( const std::vector<Identifier> & id, const std::vector<double> & hitx, const std::vector<double> & hity, const std::vector<double> & hitz, const std::vector<double> & error, const std::vector<double> & pull, std::vector<int> & select, const int n, std::vector<double> & clusterX , std::vector<double> & clusterY ,std::vector<double> & clusterZ , std::vector<double> & clusterError , std::vector<Identifier> & clusterId, std::vector<int> & clusterHits, std::vector<int> & clusterSelect, std::vector<int> & clusterInt, int & ncl ) const;
};
inline double MuonPhiHitSelector::getPhi()const{ return m_phi; }
#endif // MuonSegmentCleaner_MuonPhiHitSelector_H
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef MuonSegmentCleaner_MuonSegmentAmbiCleaner_H
#define MuonSegmentCleaner_MuonSegmentAmbiCleaner_H
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgTool.h"
#include "MuonRecToolInterfaces/IMuonSegmentCleaner.h"
class StoreGate;
class RpcIdHelper;
class MdtIdHelper;
class CscIdHelper;
class TgcIdHelper;
class StoreGateSvc;
class Identifier;
namespace Muon {
class MuonSegment;
}
namespace MuonGM {
class MuonDetectorManager;
}
namespace Trk {
class RIO_OnTrack;
class PrepRawData;
}
class MuonSegmentAmbiCleaner : public AlgTool, virtual public Muon::IMuonSegmentCleaner
{
public:
/** constructor */
MuonSegmentAmbiCleaner(const std::string&,const std::string&,const IInterface*);
/** destructor */
virtual ~MuonSegmentAmbiCleaner();
/** to initiate private members */
virtual StatusCode initialize();
/** to delete private members */
virtual StatusCode finalize();
/** For one segment solve ambiguous RPC and TGC hits: different eta but same phi
using the MDT extrapolated segment
Makes and output a new segment dropping the ambiguous hits */
virtual const Muon::MuonSegment* resolve( const Muon::MuonSegment* segment ) const;
private:
StoreGateSvc* m_storeGate; //!< Pointer to store gate
const MuonGM::MuonDetectorManager* m_detMgr; //!< Pointer to the detector manager
const RpcIdHelper* m_rpcIdHelper; //!< Pointer to RPC id helper
const TgcIdHelper* m_tgcIdHelper; //!< Pointer to TGC id helper
const CscIdHelper* m_cscIdHelper; //!< Pointer to CSC id helper
const MdtIdHelper* m_mdtIdHelper; //!< Pointer to MDT id helper
/** flag to print out debugging information */
bool m_debug;
/** flag to print out a summary of what comes in and what comes out */
bool m_summary;
};
#endif // MuonSegmentCleaner_MuonSegmentAmbiCleaner_H
package MuonSegmentCleaner
author Peter Kluit <s01@nikhef.nl>, Niels van Eldik <Niels.van.Eldik@cern.ch>
use AtlasPolicy AtlasPolicy-*
use GaudiInterface GaudiInterface-* External
use MuonRecToolInterfaces MuonRecToolInterfaces-* MuonSpectrometer/MuonReconstruction/MuonRecTools
private
use EventPrimitives EventPrimitives-* Event
use AtlasEigen AtlasEigen-* External
use CxxUtils CxxUtils-* Control
use StoreGate StoreGate-* Control
use MuonReadoutGeometry MuonReadoutGeometry-* MuonSpectrometer/MuonDetDescr
use MuonIdHelpers MuonIdHelpers-* MuonSpectrometer
use MuonPrepRawData MuonPrepRawData-* MuonSpectrometer/MuonReconstruction/MuonRecEvent
use MuonRIO_OnTrack MuonRIO_OnTrack-* MuonSpectrometer/MuonReconstruction/MuonRecEvent
use MuonSegment MuonSegment-* MuonSpectrometer/MuonReconstruction/MuonRecEvent
use MuonCompetingRIOsOnTrack MuonCompetingRIOsOnTrack-* MuonSpectrometer/MuonReconstruction/MuonRecEvent
use TrkCompetingRIOsOnTrack TrkCompetingRIOsOnTrack-* Tracking/TrkEvent
use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
use TrkPrepRawData TrkPrepRawData-* Tracking/TrkEvent
use TrkRIO_OnTrack TrkRIO_OnTrack-* Tracking/TrkEvent
use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
public
#library MuonSegmentCleaner *.cxx -s=components/*.cxx
branches MuonSegmentCleaner src
#apply_pattern component_library
apply_pattern dual_use_library files= *.cxx
#private
#macro cppdebugflags '$(cppdebugflags_s)'
#macro_remove componentshr_linkopts "-Wl,-s"
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
/**
@mainpage MuonSegmentCleaner Package
@author Peter Kluit <s01@nikhef.nl>
@author Jochem Snuverink <jsnuverink@nikhef.nl>
@section MuonSegmentCleanerIntro Introduction
This package contains tools that clean segments, e.g. outliers, delta electrons, trigger ambiguities
@section MuonSegmentCleanerOverview Class Overview
The MuonSegmentCleaner package contains the following classes:
- MuonPhiHitSelector: cleans Phi segments from outliers, build rio's, competing rio's (ambiguous hits) and perform clusterization
- MuonSegmentAmbiCleaner: solves TGC and RPC ambiguities, extrapolating MDT segments
@ref used_MuonSegmentCleaner
@ref requirements_MuonSegmentCleaner
*/
/**
@page used_MuonSegmentCleaner Used Packages
@htmlinclude used_packages.html
*/
/**
@page requirements_MuonSegmentCleaner Requirements
@include requirements
*/
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "MuonSegmentCleaner/MuonSegmentAmbiCleaner.h"
#include "MuonSegment/MuonSegment.h"
#include <sstream>
#include <iostream>
#include <vector>
#include "MuonRecToolInterfaces/IMuonClusterOnTrackCreator.h"
#include "MuonReadoutGeometry/MdtReadoutElement.h"
#include "MuonReadoutGeometry/RpcReadoutElement.h"
#include "MuonReadoutGeometry/TgcReadoutElement.h"
#include "MuonReadoutGeometry/CscReadoutElement.h"
#include "MuonReadoutGeometry/MuonDetectorManager.h"
#include "MuonIdHelpers/MdtIdHelper.h"
#include "MuonIdHelpers/RpcIdHelper.h"
#include "MuonIdHelpers/CscIdHelper.h"
#include "MuonIdHelpers/TgcIdHelper.h"
#include "MuonRIO_OnTrack/MdtDriftCircleOnTrack.h"
#include "MuonRIO_OnTrack/RpcClusterOnTrack.h"
#include "MuonRIO_OnTrack/TgcClusterOnTrack.h"
#include "MuonRIO_OnTrack/CscClusterOnTrack.h"
#include "TrkEventPrimitives/FitQuality.h"
#include "TrkEventPrimitives/LocalDirection.h"
#include "EventPrimitives/EventPrimitives.h"
#include "TrkSurfaces/PlaneSurface.h"
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "TrkPrepRawData/PrepRawData.h"
#include "TrkCompetingRIOsOnTrack/CompetingRIOsOnTrack.h"
#include "GaudiKernel/MsgStream.h"
#include "StoreGate/StoreGateSvc.h"
MuonSegmentAmbiCleaner::MuonSegmentAmbiCleaner(const std::string& type,const std::string& name,const IInterface* parent):AlgTool(type,name,parent)
{
declareInterface<IMuonSegmentCleaner>(this);
m_debug = false;
declareProperty("DoDebug",m_debug);
m_summary = false;
declareProperty("DoSummary",m_summary);
}
MuonSegmentAmbiCleaner::~MuonSegmentAmbiCleaner()
{
}
StatusCode MuonSegmentAmbiCleaner::initialize()
{
MsgStream log(msgSvc(),name());
log << MSG::VERBOSE << " MuonSegmentiAmbiCleaner::Initializing " << endreq;
StatusCode sc = service("StoreGateSvc", m_storeGate);
if (sc.isFailure()) {
log << MSG::FATAL << " StoreGate service not found " << endreq;
return StatusCode::FAILURE;
}
StoreGateSvc* detStore=0;
sc = serviceLocator()->service("DetectorStore", detStore);
if ( sc.isSuccess() ) {
sc = detStore->retrieve( m_detMgr );
if ( sc.isFailure() ) {
log << MSG::ERROR << " Cannot retrieve MuonDetDescrMgr " << endreq;
} else {
m_mdtIdHelper = m_detMgr->mdtIdHelper();
m_cscIdHelper = m_detMgr->cscIdHelper();
m_rpcIdHelper = m_detMgr->rpcIdHelper();
m_tgcIdHelper = m_detMgr->tgcIdHelper();
log << MSG::INFO << " Retrieved IdHelpers: (mdt, csc, rpc and tgc) " << endreq;
}
} else {
log << MSG::ERROR << " MuonDetDescrMgr not found in DetectorStore " << endreq;
}
log << MSG::VERBOSE << "End of Initializing" << endreq;
return StatusCode::SUCCESS;
}
StatusCode MuonSegmentAmbiCleaner::finalize()
{
return StatusCode::SUCCESS;
}
const Muon::MuonSegment* MuonSegmentAmbiCleaner::resolve(const Muon::MuonSegment* segment)const
{
MsgStream log(msgSvc(),name());
log << MSG::VERBOSE << " Executing MuonSegmentAmbiCleanerTools " << endreq;
// unsigned int nRots = segment->numberOfContainedROTs();
DataVector<const Trk::MeasurementBase>* meas_keep = new DataVector<const Trk::MeasurementBase>();
// create new surface
Trk::PlaneSurface* psf = (segment->associatedSurface()).clone();
Amg::Transform3D globalToLocal = psf->transform().inverse();
Amg::Vector3D lSegmentPos = globalToLocal*(segment->globalPosition());
Amg::Vector3D lSegmentDir = globalToLocal*(segment->globalDirection());
const std::vector<const Trk::MeasurementBase*>& meas = segment->containedMeasurements();
std::vector<const Trk::MeasurementBase*>::const_iterator mit = meas.begin();
std::vector<const Trk::MeasurementBase*>::const_iterator mit_end = meas.end();
// loop over hits
int nphirpc = 0;
int nphitgc = 0;
int nphicsc = 0;
int netamdt = 0;
int netarpc = 0;
int netatgc = 0;
int netacsc = 0;
int irio = 0;
std::vector<const Trk::RIO_OnTrack*> rots;
std::vector<const Trk::CompetingRIOsOnTrack*> crots; // lookup vector to check if rot is part of competing rio. vector contains 0 if not part of competing rio
rots.reserve(2*meas.size()); // factor 2, to be on safe side
crots.reserve(2*meas.size());
for( ; mit!=mit_end;++mit ){
//dynamic cast to either rio or competingrio:
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(*mit);
if (rot)
{
rots.push_back(rot);
crots.push_back(0);
}
else
{
const Trk::CompetingRIOsOnTrack* crio = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(*mit);
if (crio)
{
for (unsigned int i = 0; i<crio->numberOfContainedROTs(); i++)
{
rot = &crio->rioOnTrack(i);
rots.push_back(rot);
crots.push_back(crio);
}
}
}
}
unsigned int nMeas = rots.size();
unsigned int nphi = 0;
// Vectors for the phi hits
std::vector <const Trk::RIO_OnTrack*> rots_phi(nMeas);
std::vector <const Trk::CompetingRIOsOnTrack*> crots_phi(nMeas);
std::vector <const Trk::MeasurementBase*> meas_phi(nMeas);
std::vector <double> dis_phi(nMeas); // distance to segment
std::vector <int> chambercode_phi(nMeas);
std::vector <int> stripcode_phi(nMeas);
std::vector <int> ok_phi(nMeas); // 0 Not selected 1 selected/
std::vector <int> det_phi(nMeas); // 1 = RPC 2 = TGC /
std::vector <Identifier> id_phi(nMeas); // 1 = RPC 2 = TGC /
if (m_debug) std::cout << " MuonSegmentAmbiCleanerTool nMeas " << nMeas << " competing rios: " << crots.size() << std::endl;
for (unsigned int i=0; i<rots.size(); i++){
const Trk::RIO_OnTrack* rot = rots[i];
const Trk::PrepRawData* prd = rot->prepRawData();
Identifier id = prd->identify();
irio++;
// if (m_debug) std::cout << " Loop over RIOs " << irio << std::endl;
// idSegmentMap[id]= segment;
if( m_mdtIdHelper->is_mdt( rot->identify() ) ){
meas_keep->push_back(rot->clone());
netamdt++;
continue;
}else if( m_rpcIdHelper->is_rpc( rot->identify() ) ){
if( m_rpcIdHelper->measuresPhi(id) != 1) {
meas_keep->push_back(rot->clone());
netarpc++;
continue ;
}
}else if( m_tgcIdHelper->is_tgc( rot->identify() ) ){
if( m_tgcIdHelper->isStrip(id) != 1 ) {
meas_keep->push_back(rot->clone());
netatgc++;
continue;
}
}else if( m_cscIdHelper->is_csc( rot->identify() ) ){
meas_keep->push_back(rot->clone());
if( m_cscIdHelper->measuresPhi(id) != 1) {
netacsc++;
} else {
nphicsc++;
}
continue;
}
//add Phi Hits
id_phi[nphi] = id;
rots_phi[nphi] = rot;
crots_phi[nphi] = crots[i];
chambercode_phi[nphi] = 0;
stripcode_phi[nphi] = 0;
ok_phi[nphi] = 0;
det_phi[nphi] = 0;
dis_phi[nphi] = 10000000;
if (m_rpcIdHelper->is_rpc( rot->identify())) {
nphirpc++;
int code = 1000000*(m_rpcIdHelper->stationName(id));
code = code + 2*((m_rpcIdHelper->doubletR(id))-1)+16*((m_rpcIdHelper->gasGap(id))-1);
chambercode_phi[nphi] = code;
stripcode_phi[nphi] = m_rpcIdHelper->strip(id);
ok_phi[nphi] = 1;
det_phi[nphi] = 1;
const Muon::RpcClusterOnTrack* rrot = dynamic_cast<const Muon::RpcClusterOnTrack*>(rot);
if( !rrot ){
log << MSG::WARNING << "This is not a RpcClusterOnTrack!!! " << endreq;
continue;
}
const Muon::RpcPrepData* rprd = rrot->prepRawData();
Amg::Vector3D gHitPos = rprd->globalPosition();
Amg::Vector3D lHitPos = globalToLocal*gHitPos;
// In the barrel local z is measured
double disRPC = lSegmentPos.z() - lHitPos.z() + lSegmentDir.z()*(lHitPos.y()-lSegmentPos.y())/lSegmentDir.y();
if (m_debug) {
// std::cout << " lsegment pos x " << lSegmentPos.x() << " y " << lSegmentPos.y() << " z " << lSegmentPos.z() << std::endl;
// std::cout << " lsegment dir x " << lSegmentDir.x() << " y " << lSegmentDir.y() << " z " << lSegmentDir.z() << std::endl;
std::cout << " ghit pos x " << gHitPos.x() << " y " << gHitPos.y() << " z " << gHitPos.z() << std::endl;
std::cout << " dis RPC " << disRPC << std::endl;
}
dis_phi[nphi] = disRPC;
} else if ( m_tgcIdHelper->is_tgc( rot->identify())) {
nphitgc++;
int code = 1000000*(m_tgcIdHelper->stationName(id))+100*(m_tgcIdHelper->stationEta(id)+10);
code = code + m_tgcIdHelper->gasGap(id);
chambercode_phi[nphi] = code;
stripcode_phi[nphi] = m_tgcIdHelper->channel(id);
ok_phi[nphi] = 1;
det_phi[nphi] = 2;
const Muon::TgcClusterOnTrack* rrot = dynamic_cast<const Muon::TgcClusterOnTrack*>(rot);
if( !rrot ){
log << MSG::WARNING << "This is not a TgcClusterOnTrack!!! " << endreq;
continue;
}
const Muon::TgcPrepData* rprd = rrot->prepRawData();
Amg::Vector3D gHitPos = rprd->globalPosition();
Amg::Vector3D lHitPos = globalToLocal*gHitPos;
// In the forward local y is measured
double disTGC = lSegmentPos.y() - lHitPos.y() + lSegmentDir.y()*(lHitPos.z()-lSegmentPos.z())/lSegmentDir.z();
if (m_debug) {
// std::cout << " lsegment pos x " << lSegmentPos.x() << " y " << lSegmentPos.y() << " z " << lSegmentPos.z() << std::endl;
// std::cout << " lsegment dir x " << lSegmentDir.x() << " y " << lSegmentDir.y() << " z " << lSegmentDir.z() << std::endl;
std::cout << " ghit pos x " << gHitPos.x() << " y " << gHitPos.y() << " z " << gHitPos.z() << std::endl;
std::cout << " dis TGC " << disTGC << std::endl;
}
dis_phi[nphi] = disTGC;