From 423fd24c7107da1b6b51a814310fc1de944f31f1 Mon Sep 17 00:00:00 2001
From: Anthony Morley <anthony.morley@cern.ch>
Date: Thu, 27 Oct 2016 15:16:36 +0200
Subject: [PATCH] Updates for release 21 (InDetBeamSpotFinder-01-02-00)
* Update to allow:
* Fits done per unit time
* Fix additional parameters
* plb fits
---
.../InDetBeamSpotFinder/CMakeLists.txt | 16 +-
.../InDetBeamSpotFinder/IInDetBeamSpotTool.h | 1 +
.../InDetBeamSpotFinder/cmt/requirements | 47 +--
.../share/BeamspotRefitVertex.py | 19 +
.../share/postInclude.addSCTonlyToAODFile.py | 19 +
.../InDetBeamSpotFinder/src/BeamSpotID.cxx | 2 +
.../InDetBeamSpotFinder/src/BeamSpotID.h | 13 +-
.../src/InDetBeamSpotFinder.cxx | 86 ++--
.../src/InDetBeamSpotFinder.h | 3 +
.../src/InDetBeamSpotReader.cxx | 47 +--
.../src/InDetBeamSpotRooFit.cxx | 12 +-
.../src/InDetBeamSpotVertex.cxx | 363 +++++++++--------
.../src/InDetBeamSpotVertex.h | 7 +-
.../src/RefitTracksAndVertex.cxx | 382 ++++++++++++++++++
.../src/RefitTracksAndVertex.h | 90 +++++
.../InDetBeamSpotFinder_entries.cxx | 5 +
16 files changed, 837 insertions(+), 275 deletions(-)
create mode 100644 InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/BeamspotRefitVertex.py
create mode 100644 InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/postInclude.addSCTonlyToAODFile.py
create mode 100644 InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.cxx
create mode 100644 InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.h
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/CMakeLists.txt b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/CMakeLists.txt
index 46dcb530d582..35df023212b0 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/CMakeLists.txt
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/CMakeLists.txt
@@ -16,7 +16,17 @@ atlas_depends_on_subdirs( PUBLIC
InnerDetector/InDetConditions/InDetBeamSpotService
Tracking/TrkEvent/TrkEventPrimitives
Tracking/TrkEvent/VxVertex
- Trigger/TrigAnalysis/TrigAnalysisInterfaces )
+ Trigger/TrigAnalysis/TrigAnalysisInterfaces
+ DetectorDescription/AtlasDetDescr
+ DetectorDescription/Identifier
+ Tracking/TrkFitter/TrkFitterInterfaces
+ Tracking/TrkFitter/TrkFitterUtils
+ Tracking/TrkEvent/TrkParameters
+ Tracking/TrkEvent/TrkPseudoMeasurementOnTrack
+ Tracking/TrkEvent/TrkRIO_OnTrack
+ Tracking/TrkEvent/TrkTrack
+ Tracking/TrkVertexFitter/TrkVertexFitterInterfaces
+)
# External dependencies:
find_package( CLHEP )
@@ -30,12 +40,12 @@ atlas_add_library( InDetBeamSpotFinderLib
PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
DEFINITIONS ${CLHEP_DEFINITIONS}
LINK_LIBRARIES ${CLHEP_LIBRARIES} xAODTracking GaudiKernel
- PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps EventInfo xAODEventInfo TrkEventPrimitives VxVertex )
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps EventInfo xAODEventInfo TrkEventPrimitives VxVertex AtlasDetDescr Identifier TrkFitterInterfaces TrkFitterUtils TrkParameters TrkPseudoMeasurementOnTrack TrkRIO_OnTrack TrkTrack TrkVertexFitterInterfaces )
atlas_add_component( InDetBeamSpotFinder
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} xAODTracking GaudiKernel AthenaBaseComps EventInfo xAODEventInfo TrkEventPrimitives VxVertex InDetBeamSpotFinderLib )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} xAODTracking GaudiKernel AthenaBaseComps EventInfo xAODEventInfo TrkEventPrimitives VxVertex AtlasDetDescr Identifier TrkFitterInterfaces TrkFitterUtils TrkParameters TrkPseudoMeasurementOnTrack TrkRIO_OnTrack TrkTrack TrkVertexFitterInterfaces InDetBeamSpotFinderLib )
# Install files from the package:
atlas_install_joboptions( share/*.py )
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/InDetBeamSpotFinder/IInDetBeamSpotTool.h b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/InDetBeamSpotFinder/IInDetBeamSpotTool.h
index 55bdee2b0a5e..7c827f51c979 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/InDetBeamSpotFinder/IInDetBeamSpotTool.h
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/InDetBeamSpotFinder/IInDetBeamSpotTool.h
@@ -31,6 +31,7 @@ namespace BeamSpot {
namespace BeamSpot{
struct Event {
unsigned int pileup, runNumber, lumiBlock, bcid;
+ unsigned long long eventNumber, eventTime, eventTime_NS;
std::vector<BeamSpot::VrtHolder> vertices;
};
}
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/cmt/requirements b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/cmt/requirements
index 0e248367f090..69fc53d5b293 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/cmt/requirements
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/cmt/requirements
@@ -4,40 +4,33 @@ author James Walder <jwalder@cern.ch>
author Guennadi Borissov <Gennady.Borisov@cern.ch>
author Brian Amadio <btamadio@lbl.gov>
-use AtlasPolicy AtlasPolicy-*
-use GaudiInterface GaudiInterface-* External
-#use StoreGate StoreGate-* Control
-#private
-#use AthenaKernel AthenaKernel-* Control
-public
-use AtlasCLHEP AtlasCLHEP-* External
+use AtlasPolicy AtlasPolicy-*
+use GaudiInterface GaudiInterface-* External
+use AtlasCLHEP AtlasCLHEP-* External
+use xAODTracking xAODTracking-* Event/xAOD
private
-use AthenaBaseComps AthenaBaseComps-* Control
-use VxVertex VxVertex-* Tracking/TrkEvent
-#use GeneratorObjects GeneratorObjects-* Generators
+use AthenaBaseComps AthenaBaseComps-* Control
+use VxVertex VxVertex-* Tracking/TrkEvent
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
-#public
-#use Particle Particle-* Reconstruction
-#use TrkTrack TrkTrack-* Tracking/TrkEvent
-#use TrkToolInterfaces TrkToolInterfaces-* Tracking/TrkTools
+use EventInfo EventInfo-* Event
+use xAODEventInfo xAODEventInfo-* Event/xAOD
-private
-use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
-#use TrkTrackSummary TrkTrackSummary-* Tracking/TrkEvent
-#use TrkParameters TrkParameters-* Tracking/TrkEvent
+use InDetBeamSpotService InDetBeamSpotService-* InnerDetector/InDetConditions
+use TrigAnalysisInterfaces TrigAnalysisInterfaces-* Trigger/TrigAnalysis
-private
-use EventInfo EventInfo-* Event
-use xAODEventInfo xAODEventInfo-* Event/xAOD
+use AtlasDetDescr AtlasDetDescr-* DetectorDescription
+use Identifier Identifier-* DetectorDescription
+use TrkFitterInterfaces TrkFitterInterfaces-* Tracking/TrkFitter
+use TrkFitterUtils TrkFitterUtils-* Tracking/TrkFitter
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkPseudoMeasurementOnTrack TrkPseudoMeasurementOnTrack-* Tracking/TrkEvent
+use TrkRIO_OnTrack TrkRIO_OnTrack-* Tracking/TrkEvent
+use TrkTrack TrkTrack-* Tracking/TrkEvent
+use TrkVertexFitterInterfaces TrkVertexFitterInterfaces-* Tracking/TrkVertexFitter
-public
-use xAODTracking xAODTracking-* Event/xAOD
-
-private
-use InDetBeamSpotService InDetBeamSpotService-* InnerDetector/InDetConditions
-use TrigAnalysisInterfaces TrigAnalysisInterfaces-* Trigger/TrigAnalysis
private
apply_tag ROOTMathLibs
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/BeamspotRefitVertex.py b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/BeamspotRefitVertex.py
new file mode 100644
index 000000000000..f197e7b206da
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/BeamspotRefitVertex.py
@@ -0,0 +1,19 @@
+
+# Example
+#Reco_tf.py --steering 'doRAWtoALL' --ignoreErrors 'True' --inputBSFile /tmp/amorley/data16_13TeV.00310809.calibration_BeamSpot.merge.RAW/data16_13TeV.00310809.calibration_BeamSpot.merge.RAW._lb0743._SFO-ALL._0001.1 --outputAODFile /tmp/amorley/AOD.pool.root --conditionsTag all:CONDBR2-ES1PA-2016-03 --AMITag 'c1042' --autoConfiguration='everything' --maxEvents '-1' --preInclude 'all:InDetBeamSpotExample/preIncludeRecoForBeamspot.py' --postExec 'all:from IOVDbSvc.CondDB import conddb; conddb.addOverride("/Indet/AlignL1/ID" ,"IndetAlignL1ID-RUN2-BLK-UPD4-02"); conddb.addOverride("/Indet/AlignL2/PIX" ,"IndetAlignL2PIX-RUN2-BLK-UPD4-02"); conddb.addOverride("/Indet/AlignL2/SCT" ,"IndetAlignL2SCT-RUN2-BLK-UPD4-02"); conddb.addOverride("/Indet/AlignL3" ,"IndetAlignL3-RUN2-BLK-UPD4-02"); conddb.addOverride("/Indet/IBLDist" ,"InDetIBLDist-RUN2-BLK-UPD4-01"); conddb.addOverride("/TRT/AlignL1/TRT" ,"TRTAlignL1-RUN2-BLK-UPD4-02"); conddb.addOverride("/TRT/AlignL2" ,"TRTAlignL2-RUN2-BLK-UPD4-02");' --preExec 'all:from InDetRecExample.InDetJobProperties import InDetFlags; InDetFlags.useDynamicAlignFolders.set_Value_and_Lock(True); from RecExConfig.RecFlags import rec;rec.UserAlgs.set_Value_and_Lock("InDetBeamSpotFinder/BeamspotRefitVertex.py")' --postInclude 'all:BeamSpotFinder/postInclude.addSCTonlyToAODFile.py' --geometryVersion all:ATLAS-R2-2015-04-00-00 --maxEvents 10
+
+from InDetRecExample.InDetKeys import InDetKeys
+containerName= "SCTonlyVertex"
+from InDetBeamSpotFinder.InDetBeamSpotFinderConf import RefitTracksAndVertex
+thisRefitTracksAndVertex = RefitTracksAndVertex(name = 'RefitTracksAndVertex',
+ TrackFitter = ToolSvc.InDetTrackFitter,
+ VertexFitterTool = InDetVxFitterTool,
+ VertexListInput = InDetKeys.xAODVertexContainer(),
+ OutputVertexContainer = containerName)
+#thisRefitTracksAndVertex.SelEtaMin = -0.8
+#thisRefitTracksAndVertex.SelEtaMax = 0.8
+#thisRefitTracksAndVertex.OutputLevel = 1
+
+from RecExConfig.RecFlags import rec
+topSequence += thisRefitTracksAndVertex
+
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/postInclude.addSCTonlyToAODFile.py b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/postInclude.addSCTonlyToAODFile.py
new file mode 100644
index 000000000000..563d1f7188e6
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/postInclude.addSCTonlyToAODFile.py
@@ -0,0 +1,19 @@
+containerName= "SCTonlyVertex"
+AnAODList = []
+AnAODList+=['xAOD::VertexContainer#'+containerName]
+AnAODList+=['xAOD::VertexAuxContainer#'+containerName+'Aux.-vxTrackAtVertex']
+
+StreamAOD.ItemList += AnAODList
+
+
+
+if False:
+ from OutputStreamAthenaPool.MultipleStreamManager import MSMgr
+ SCTonlyVertexStream = MSMgr.NewPoolRootStream( "SCTonlyVertexStream", "SCTonlyVertex.pool.root" )
+ SCTonlyVertexStream.AddItem( "TrigDec::TrigDecision#TrigDecision" )
+ SCTonlyVertexStream.AddItem( 'xAOD::VertexContainer#'+containerName )
+ SCTonlyVertexStream.AddItem( 'xAOD::VertexAuxContainer#'+containerName+'Aux.-vxTrackAtVertex' )
+ SCTonlyVertexStream.AddItem( 'xAOD::EventInfo#EventInfo' )
+ SCTonlyVertexStream.AddItem( 'xAOD::EventAuxInfo#EventInfoAux.' )
+
+
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.cxx
index 6e054897994b..2b4e7c94cc8d 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.cxx
@@ -5,6 +5,8 @@
#include "BeamSpotID.h"
bool BeamSpot::ID::operator<( const ID & rhs) const {
+ if ( m_timeStamp < rhs.m_timeStamp ) return true;
+ if ( m_timeStamp > rhs.m_timeStamp ) return false;
if ( m_pileup < rhs.m_pileup ) return true;
if ( m_pileup > rhs.m_pileup ) return false;
if ( m_bcid < rhs.m_bcid ) return true;
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.h b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.h
index bbb5e57f6796..6d9041ee88f2 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.h
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/BeamSpotID.h
@@ -17,20 +17,23 @@ namespace BeamSpot {
unsigned int lumiBlock() const {return m_lumiBlock; }
unsigned int pileup() const {return m_pileup; }
unsigned int bcid() const {return m_bcid; }
+ unsigned long timeStamp() const{return m_timeStamp; }
void runNumber(unsigned int run) {m_runNumber = run;}
void lumiBlock(unsigned int lb ) {m_lumiBlock = lb;}
void pileup(unsigned int pileup) {m_pileup = pileup;}
void bcid(unsigned int bcid) {m_bcid = bcid; }
-
+ void timeStamp(unsigned long time){m_timeStamp = time;}
+
bool operator<( const ID & ) const;
private:
- unsigned int m_runNumber;
- unsigned int m_lumiBlock;
- unsigned int m_pileup;
- unsigned int m_bcid;
+ unsigned int m_runNumber;
+ unsigned int m_lumiBlock;
+ unsigned int m_pileup;
+ unsigned int m_bcid;
+ unsigned long m_timeStamp;
};
}
#endif
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.cxx
index a16296f79fd8..47da56076488 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.cxx
@@ -12,6 +12,9 @@
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/TrackParticle.h"
+#include "EventInfo/EventInfo.h"
+#include "EventInfo/EventID.h"
+
InDet::InDetBeamSpotFinder::InDetBeamSpotFinder(const std::string& name, ISvcLocator* pSvcLocator):
AthAlgorithm(name, pSvcLocator),
m_toolSvc("ToolSvc",name),m_bcTool("Trig::TrigConfBunchCrossingTool/BunchCrossingTool")
@@ -36,12 +39,13 @@ InDet::InDetBeamSpotFinder::InDetBeamSpotFinder(const std::string& name, ISvcLoc
declareProperty( "VertexNtuple" , m_writeVertexNtuple = true);
declareProperty( "WriteAllVertices" , m_writeAllVertices=false);
declareProperty( "VertexTreeName" , m_vertexTreeName);
+ declareProperty( "SecondsPerFit", m_secondsPerFit = 1);
}
StatusCode InDet::InDetBeamSpotFinder::initialize() {
- msg(MSG::DEBUG) << "in initialize()" << endreq;
+ ATH_MSG_DEBUG( "in initialize()");
if ( m_beamSpotToolList.size() == 0 ){
- msg(MSG::FATAL) <<"FATAL ERROR: must provide at least one beamspot tool in beamSpotToolList";
+ ATH_MSG_FATAL("FATAL ERROR: must provide at least one beamspot tool in beamSpotToolList");
return StatusCode::FAILURE;
}
ATH_CHECK( service("THistSvc",m_thistSvc) );
@@ -65,7 +69,7 @@ StatusCode InDet::InDetBeamSpotFinder::execute(){
if( m_writeVertexNtuple ){
for( unsigned int i = 0; i < currentEvent.vertices.size(); i++){
if( currentEvent.vertices[i].passed || m_writeAllVertices ){
- writeToVertexTree( currentEvent, currentEvent.vertices[i] );
+ writeToVertexTree( currentEvent, currentEvent.vertices[i] );
}
}
}
@@ -73,7 +77,7 @@ StatusCode InDet::InDetBeamSpotFinder::execute(){
}
StatusCode InDet::InDetBeamSpotFinder::finalize() {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "in finalize()" << endreq;
+ ATH_MSG_DEBUG( "in finalize()" );
sortEvents();
ATH_CHECK(performFits());
return StatusCode::SUCCESS;
@@ -85,7 +89,22 @@ BeamSpot::Event InDet::InDetBeamSpotFinder::readEvent(const xAOD::EventInfo & ev
event.pileup = ceil(eventInfo.actualInteractionsPerCrossing());
event.runNumber = eventInfo.runNumber();
event.lumiBlock = eventInfo.lumiBlock();
- event.bcid = eventInfo.bcid();
+ event.bcid = eventInfo.bcid();
+ event.eventTime = eventInfo.timeStamp();
+ event.eventTime_NS = eventInfo.timeStampNSOffset();
+ event.eventNumber = eventInfo.eventNumber();
+ const DataHandle<EventInfo> BSeventInfo;
+ //This is required for pseudo lumiblocks
+ if( evtStore()->retrieve(BSeventInfo) != StatusCode::SUCCESS){
+ ATH_MSG_ERROR("Cannot get event info.");
+ return event;
+ }
+ if (event.lumiBlock != BSeventInfo->event_ID()->lumi_block())
+ {
+ //ATH_MSG_INFO("Updating Event info " << BSeventInfo->event_ID()->lumi_block() );
+ event.lumiBlock = BSeventInfo->event_ID()->lumi_block();
+ }
+
for(xAOD::VertexContainer::const_iterator vtxItr=vertexContainer.begin(); vtxItr!=vertexContainer.end(); ++vtxItr) {
if ((*vtxItr)->vertexType() == xAOD::VxType::NoVtx) continue;
vertex.x = (*vtxItr)->x();
@@ -111,6 +130,8 @@ void InDet::InDetBeamSpotFinder::sortEvents(){
id.lumiBlock( (m_maxLBsPerFit > 0) ? m_eventList[i].lumiBlock : 0 );
id.pileup ( iequals(m_fitSortingKey,"pileup") ? m_eventList[i].pileup : 0 );
id.bcid ( iequals(m_fitSortingKey,"bcid" ) ? m_eventList[i].bcid : 0 );
+ //std::cout << "time " << iequals(m_fitSortingKey,"time" )<< " " << m_eventList[i].eventTime/m_secondsPerFit << std::endl;
+ id.timeStamp( iequals(m_fitSortingKey,"time" ) ? m_eventList[i].eventTime/m_secondsPerFit : 0 );
m_eventMap[id].push_back( m_eventList[i] );
}
auto iter = m_eventMap.begin();
@@ -124,18 +145,19 @@ void InDet::InDetBeamSpotFinder::sortEvents(){
currentID = iter->first;
if( iter == m_eventMap.begin() || currentID.runNumber() != lastID.runNumber() ){ nRuns++; }
if( iter == m_eventMap.begin() || currentID.lumiBlock() != lastID.lumiBlock() ){ nLBs++; }
- if( currentID.pileup() != lastID.pileup() || currentID.bcid() != lastID.bcid()
- || ( m_maxRunsPerFit > 0 && nRuns > m_maxRunsPerFit )
- || ( m_maxLBsPerFit > 0 && nLBs > m_maxLBsPerFit )){
+ if( currentID.timeStamp() != lastID.timeStamp() ||
+ currentID.pileup() != lastID.pileup() || currentID.bcid() != lastID.bcid()
+ || ( m_maxRunsPerFit > 0 && nRuns > m_maxRunsPerFit )
+ || ( m_maxLBsPerFit > 0 && nLBs > m_maxLBsPerFit )){
nFits++;
m_sortedEventList.resize(nFits);
nRuns = 1; nLBs = 1;
}
for( unsigned int i = 0; i < iter->second.size(); i++){
if( m_sortedEventList.at(nFits-1).size() == m_maxEventsPerFit && m_maxEventsPerFit > 0 ){
- nFits++;
- m_sortedEventList.resize(nFits);
- nRuns = 1; nLBs = 1;
+ nFits++;
+ m_sortedEventList.resize(nFits);
+ nRuns = 1; nLBs = 1;
}
m_sortedEventList.at(nFits-1).push_back( iter->second.at(i) );
}
@@ -147,7 +169,7 @@ void InDet::InDetBeamSpotFinder::sortEvents(){
void InDet::InDetBeamSpotFinder::convertVtxTypeNames(){
if ( m_vertexTypeNames.size() ) {
for ( std::vector<std::string>::const_iterator it = m_vertexTypeNames.begin();
- it != m_vertexTypeNames.end(); ++it) {
+ it != m_vertexTypeNames.end(); ++it) {
if ((*it) == "NoVtx") ;
else if ((*it) == "PriVtx") m_vertexTypes.push_back(xAOD::VxType::PriVtx);
else if ((*it) == "SecVtx") m_vertexTypes.push_back(xAOD::VxType::SecVtx);
@@ -157,11 +179,10 @@ void InDet::InDetBeamSpotFinder::convertVtxTypeNames(){
else if ((*it) == "KinkVtx") m_vertexTypes.push_back(xAOD::VxType::KinkVtx);
else if ((*it) =="NotSpecified")m_vertexTypes.push_back(xAOD::VxType::NotSpecified) ;
}
- msg(MSG::INFO) << "Allowing " << m_vertexTypes.size()
- << " Vertex types" << endreq;
+ ATH_MSG_INFO("Allowing " << m_vertexTypes.size() << " Vertex types" );
}
else {
- msg(MSG::DEBUG) << "No selection based on vertexType will be done" << endreq;
+ ATH_MSG_DEBUG( "No selection based on vertexType will be done" );
}
}
@@ -192,7 +213,7 @@ StatusCode InDet::InDetBeamSpotFinder::setupVertexTree(){
const std::string inRootID = "/INDETBEAMSPOTFINDER/";
const std::string svrts = m_vertexTreeName;
m_root_vrt = new TTree(svrts.data(),"Vertices");
- m_root_vrt->Branch("vrt",&m_root_vtx,"x/D:y:z:vxx:vxy:vyy:vzz:vType/i:run:lb:bcid:pileup:nTracks:passed/O:valid");
+ m_root_vrt->Branch("vrt",&m_root_vtx,"x/D:y:z:vxx:vxy:vyy:vzz:vType/i:run:lb:bcid:pileup:nTracks:eventNumber/l:eventTime:eventTime_NS:passed/O:valid");
return m_thistSvc->regTree(inRootID+svrts,m_root_vrt);
}
@@ -204,7 +225,7 @@ StatusCode InDet::InDetBeamSpotFinder::performFits(){
verticesToFit.clear();
for( unsigned int j = 0; j < m_sortedEventList.at(i).size(); j++){
for( unsigned int k = 0; k < m_sortedEventList.at(i).at(j).vertices.size(); k++){
- if( m_sortedEventList.at(i).at(j).vertices.at(k).passed ) { verticesToFit.push_back( m_sortedEventList.at(i).at(j).vertices.at(k) ); }
+ if( m_sortedEventList.at(i).at(j).vertices.at(k).passed ) { verticesToFit.push_back( m_sortedEventList.at(i).at(j).vertices.at(k) ); }
}
}
for( unsigned int j = 0; j < m_beamSpotToolList.size(); j++){
@@ -219,19 +240,19 @@ StatusCode InDet::InDetBeamSpotFinder::performFits(){
m_BeamStatusCode.setAlgType( bs->getFitID() );
int fitStat;
if ( bsFitStatus == IInDetBeamSpotTool::successful)
- fitStat = 3;
+ fitStat = 3;
else if ( bsFitStatus == IInDetBeamSpotTool::problems)
- fitStat = 1;
+ fitStat = 1;
else if ( bsFitStatus == IInDetBeamSpotTool::failed || bsFitStatus == IInDetBeamSpotTool::unsolved)
- fitStat = 0;
+ fitStat = 0;
else
- fitStat = 0;
+ fitStat = 0;
m_BeamStatusCode.setFitStatus(fitStat);
if (bs->getParamMap()["sigmaX"] == 0 && bs->getParamMap()["sigmaY"] ==0 ) { m_BeamStatusCode.setFitWidth( false); }
else { m_BeamStatusCode.setFitWidth(true); }
if(m_sortedEventList[i].size() > 0)
- writeToBeamSpotTree( bs, m_sortedEventList[i], verticesToFit );
+ writeToBeamSpotTree( bs, m_sortedEventList[i], verticesToFit );
}
}
return StatusCode::SUCCESS;
@@ -270,9 +291,9 @@ StatusCode InDet::InDetBeamSpotFinder::setupBeamSpotTree(){
std::string key = iter->first;
//double val = iter->second;
if( !(m_root_bs->GetBranch(key.c_str())) ){
- m_beamSpotNtuple.paramMap[key] = 0;
- keySlashD = key + slashD;
- m_root_bs->Branch(key.c_str(), &m_beamSpotNtuple.paramMap[key], keySlashD.c_str());
+ m_beamSpotNtuple.paramMap[key] = 0;
+ keySlashD = key + slashD;
+ m_root_bs->Branch(key.c_str(), &m_beamSpotNtuple.paramMap[key], keySlashD.c_str());
}
}
//Loop over the covariance matrix for a given fit tool and create a branch for each element, if it doesn't already exist.
@@ -280,9 +301,9 @@ StatusCode InDet::InDetBeamSpotFinder::setupBeamSpotTree(){
std::string key = iter->first;
//double val = iter->second;
if( !(m_root_bs->GetBranch(key.c_str())) ){
- m_beamSpotNtuple.covMap[key] = 0;
- keySlashD = key + slashD;
- m_root_bs->Branch( key.c_str(), &m_beamSpotNtuple.covMap[key], keySlashD.c_str());
+ m_beamSpotNtuple.covMap[key] = 0;
+ keySlashD = key + slashD;
+ m_root_bs->Branch( key.c_str(), &m_beamSpotNtuple.covMap[key], keySlashD.c_str());
}
}
}
@@ -304,7 +325,10 @@ void InDet::InDetBeamSpotFinder::writeToVertexTree(BeamSpot::Event &evt, BeamSpo
m_root_vtx.pileup = evt.pileup;
m_root_vtx.nTracks = vtx.nTracks;
m_root_vtx.passed = vtx.passed;
- m_root_vtx.valid = vtx.valid;
+ m_root_vtx.valid = vtx.valid;
+ m_root_vtx.eventNumber = evt.eventNumber;
+ m_root_vtx.eventTime = evt.eventTime;
+ m_root_vtx.eventTime_NS = evt.eventTime_NS;
m_root_vrt->Fill();
}
@@ -341,8 +365,8 @@ void InDet::InDetBeamSpotFinder::writeToBeamSpotTree(const IInDetBeamSpotTool *b
m_beamSpotNtuple.run = min_run( eventList );
m_beamSpotNtuple.status = m_BeamStatusCode.getWord();
m_beamSpotNtuple.separation = 0;
- m_beamSpotNtuple.timeEnd = 0;
- m_beamSpotNtuple.timeStart = 0;
+ m_beamSpotNtuple.timeEnd = iequals(m_fitSortingKey,"time") ? eventList.back().eventTime : 0;
+ m_beamSpotNtuple.timeStart = iequals(m_fitSortingKey,"time") ? eventList.front().eventTime : 0;
m_beamSpotNtuple.runEnd = max_run( eventList );
for( std::map<std::string,double>::iterator iter = m_beamSpotNtuple.paramMap.begin();
iter != m_beamSpotNtuple.paramMap.end(); ++iter){
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.h b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.h
index a87e48e02cfb..0af945fed78c 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.h
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotFinder.h
@@ -68,6 +68,7 @@ namespace InDet {
double x,y,z,vxx,vxy,vyy,vzz;
xAOD::VxType::VertexType vType;
unsigned int run, lb, bcid, pileup, nTracks;
+ unsigned long long eventNumber, eventTime, eventTime_NS;
bool passed, valid;
};
@@ -134,6 +135,8 @@ namespace InDet {
//std::map<BeamSpot::ID, long> m_nEvents;
unsigned int m_pileupMin;
unsigned int m_pileupMax;
+
+ unsigned long m_secondsPerFit;
};
}//end namespace
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotReader.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotReader.cxx
index 1596be7cb94e..ce46a8bdcfaa 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotReader.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotReader.cxx
@@ -29,21 +29,21 @@ InDet::InDetBeamSpotReader::InDetBeamSpotReader(const std::string& name, ISvcLoc
StatusCode InDet::InDetBeamSpotReader::initialize() {
- msg(MSG::DEBUG) << "in initialize()" << endreq;
+ ATH_MSG_DEBUG( "in initialize()" );
if ( m_toolSvc.retrieve().isFailure() ) {
- msg(MSG::FATAL) << "Failed to retrieve service " << m_toolSvc << endreq;
+ ATH_MSG_FATAL( "Failed to retrieve service " << m_toolSvc );
return StatusCode::FAILURE;
} else
- msg(MSG::INFO) << "Retrieved service " << m_toolSvc << endreq;
+ ATH_MSG_INFO( "Retrieved service " << m_toolSvc );
if ( m_beamSpotSvc.retrieve().isFailure() ) {
- msg(MSG::FATAL) << "Failed to retrieve service " << m_beamSpotSvc << endreq;
+ ATH_MSG_FATAL( "Failed to retrieve service " << m_beamSpotSvc );
return StatusCode::FAILURE;
} else
- msg(MSG::INFO) << "Retrieved service " << m_beamSpotSvc << endreq;
+ ATH_MSG_INFO( "Retrieved service " << m_beamSpotSvc );
return StatusCode::SUCCESS;
@@ -51,36 +51,31 @@ StatusCode InDet::InDetBeamSpotReader::initialize() {
StatusCode InDet::InDetBeamSpotReader::execute(){
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "in execute()" << endreq;
+ ATH_MSG_DEBUG( "in execute()");
//get the set of
const DataHandle<EventInfo> eventInfo;
if (StatusCode::SUCCESS != evtStore()->retrieve( eventInfo ) ){
- if (msgLvl(MSG::ERROR)) msg(MSG::ERROR) << "Cannot get event info." << endreq;
+ ATH_MSG_ERROR( "Cannot get event info." );
return StatusCode::FAILURE;
}
EventID* eventID = eventInfo->event_ID();
- if (msgLvl(MSG::INFO)) {
- msg(MSG::INFO) << "In event " << (*eventID) << endreq;
-
- msg(MSG::INFO) <<"BeamSpot Position: \n "
- << m_beamSpotSvc->beamPos() << endreq;
- msg(MSG::INFO) <<"BeamSpot Sigma\n\t"
+ ATH_MSG_INFO( "In event " << (*eventID) );
+ ATH_MSG_INFO("BeamSpot Position: \n "
+ << m_beamSpotSvc->beamPos() );
+ ATH_MSG_INFO("BeamSpot Sigma\n\t"
<< m_beamSpotSvc->beamSigma(0) << "\n\t"
<< m_beamSpotSvc->beamSigma(1) << "\n\t"
- << m_beamSpotSvc->beamSigma(2) << "\n\t"
- << endreq;
- msg(MSG::INFO) <<"BeamSpot Tilt\n\t"
+ << m_beamSpotSvc->beamSigma(2) << "\n\t");
+ ATH_MSG_INFO("BeamSpot Tilt\n\t"
<< m_beamSpotSvc->beamTilt(0) << "\n\t"
- << m_beamSpotSvc->beamTilt(1) << "\n\t"
- << endreq;
- msg(MSG::INFO) <<"Beamspot position at PV z-position" << endreq;
- }
+ << m_beamSpotSvc->beamTilt(1) << "\n\t");
+ ATH_MSG_INFO("Beamspot position at PV z-position");
const VxContainer* importedVxContainer =0;
static const std::string m_containerName = "VxPrimaryCandidate";
if ( StatusCode::SUCCESS != evtStore()->retrieve(importedVxContainer,m_containerName)){
- if (msgLvl(MSG::ERROR)) msg(MSG::ERROR) << "No " << m_containerName << " found in StoreGate" << endreq;
+ ATH_MSG_ERROR( "No " << m_containerName << " found in StoreGate" );
return StatusCode::FAILURE;
}
//get list of PVs
@@ -88,16 +83,16 @@ StatusCode InDet::InDetBeamSpotReader::execute(){
for(vtxItr=importedVxContainer->begin();
vtxItr!=importedVxContainer->end(); ++vtxItr) {
if (static_cast<int>((*vtxItr)->vxTrackAtVertex()->size())==0) continue;
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) <<"PV position: "
- << (*vtxItr)->recVertex().position()<< endreq;
+ if (msgLvl(MSG::INFO)) ATH_MSG_INFO("PV position: "
+ << (*vtxItr)->recVertex().position() );
double z = (*vtxItr)->recVertex().position().z();
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) <<"\n\t"
+ if (msgLvl(MSG::INFO)) ATH_MSG_INFO("\n\t"
<< m_beamSpotSvc->beamPos()(0)
+ (z - m_beamSpotSvc->beamPos()(2))
*m_beamSpotSvc->beamTilt(0) << "\n\t"
<< m_beamSpotSvc->beamPos()(1)
+ (z - m_beamSpotSvc->beamPos()(2))
- *m_beamSpotSvc->beamTilt(1) << endreq;
+ *m_beamSpotSvc->beamTilt(1) );
}
return StatusCode::SUCCESS;
@@ -106,7 +101,7 @@ StatusCode InDet::InDetBeamSpotReader::execute(){
StatusCode InDet::InDetBeamSpotReader::finalize() {
- msg(MSG::DEBUG) << "in finalize()" << endreq;
+ ATH_MSG_DEBUG( "in finalize()" );
return StatusCode::SUCCESS;
}
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotRooFit.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotRooFit.cxx
index 17b479615ea5..dc21a4702438 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotRooFit.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotRooFit.cxx
@@ -48,13 +48,12 @@ InDetBeamSpotRooFit::InDetBeamSpotRooFit( const InDetBeamSpotRooFit& rhs ) :
{}
StatusCode InDetBeamSpotRooFit::initialize() {
- msg().setLevel( outputLevel() );
- msg(MSG::DEBUG) << "In initialize" << endreq;
+ ATH_MSG_DEBUG( "In initialize" );
return StatusCode::SUCCESS;
}
StatusCode InDetBeamSpotRooFit::finalize() {
- msg(MSG::DEBUG) << "In finalize" << endreq;
+ ATH_MSG_DEBUG( "In finalize");
return StatusCode::SUCCESS;
}
@@ -171,7 +170,7 @@ InDetBeamSpotRooFit::FitStatus InDetBeamSpotRooFit::fit(std::vector< BeamSpot::V
std::string combinedCutString = m_vtxCutString + std::string(" && ") + rmsCutString;
const char *combinedCut = (const char*)combinedCutString.c_str();
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) <<"combinedCut = "<<combinedCut<<endreq;
+ ATH_MSG_INFO( "combinedCut = "<<combinedCut );
//Convert these numbers to strings, then concatenate the strings
@@ -184,9 +183,8 @@ InDetBeamSpotRooFit::FitStatus InDetBeamSpotRooFit::fit(std::vector< BeamSpot::V
r->Print();
m_nUsed = rfData.reduce(Cut(combinedCut))->numEntries();
- if (msgLvl(MSG::INFO)) msg(MSG::INFO)
- << "A total of " << m_vertexCount << " vertices passed pre-selection. Of which "
- << rfData.reduce(Cut(combinedCut))->numEntries()<<" vertices will be used in the ML fit."<<endreq;
+ ATH_MSG_INFO( "A total of " << m_vertexCount << " vertices passed pre-selection. Of which "
+ << rfData.reduce(Cut(combinedCut))->numEntries()<<" vertices will be used in the ML fit.") ;
if ( r->edm() <= 10e-04 && r->covQual() == 3){
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.cxx
index 370320a1d40f..67c1f76b47cb 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.cxx
@@ -36,8 +36,8 @@ public:
InDetBeamSpotVertex::InDetBeamSpotVertex( const std::string& type,
- const std::string& name,
- const IInterface* parent):
+ const std::string& name,
+ const IInterface* parent):
AthAlgTool(type,name,parent),
m_x(4),m_cov(4,0),m_z(0.),m_zErr(0.),
m_p(4), m_V(4,0),
@@ -56,8 +56,8 @@ InDetBeamSpotVertex::InDetBeamSpotVertex( const std::string& type,
declareProperty("InitParZ", m_init_z =0.0);
declareProperty("InitParAx", m_init_ax=0.0);
declareProperty("InitParAy", m_init_ay=0.0);
- declareProperty("InitParSigmaX", m_init_sx=0.015);
- declareProperty("InitParSigmaY", m_init_sy=0.015);
+ declareProperty("InitParSigmaX", m_init_sx=0.01);
+ declareProperty("InitParSigmaY", m_init_sy=0.01);
declareProperty("InitParSigmaZ", m_init_sz=56.0);
declareProperty("InitParK", m_init_k =1.0);
declareProperty("InitParRhoXY", m_init_rhoxy=0.0);
@@ -67,8 +67,8 @@ InDetBeamSpotVertex::InDetBeamSpotVertex( const std::string& type,
declareProperty("InitParMinZ", m_init_min_z =0.0);
declareProperty("InitParMinAx", m_init_min_ax=0.0);
declareProperty("InitParMinAy", m_init_min_ay=0.0);
- declareProperty("InitParMinSigmaX", m_init_min_sx=0.001);
- declareProperty("InitParMinSigmaY", m_init_min_sy=0.001);
+ declareProperty("InitParMinSigmaX", m_init_min_sx=0.0001);
+ declareProperty("InitParMinSigmaY", m_init_min_sy=0.0001);
declareProperty("InitParMinSigmaZ", m_init_min_sz=0.0);
declareProperty("InitParMinK", m_init_min_k =0);
@@ -114,6 +114,7 @@ InDetBeamSpotVertex::InDetBeamSpotVertex( const std::string& type,
declareProperty( "FixParK" , m_fixInputK = false);
declareProperty( "UseLLNorm" , m_useLLNorm = false);
+ declareProperty( "FixWidth", m_fixWidth = false);
declareProperty("TruncatedRMS", m_truncatedRMS = false);
declareProperty("RMSFraction", m_fractionRMS = 0.95);
@@ -197,6 +198,7 @@ InDetBeamSpotVertex::InDetBeamSpotVertex( const InDetBeamSpotVertex& rhs) :
m_fixInputK = rhs.m_fixInputK;
m_useLLNorm = rhs.m_useLLNorm;
+ m_fixWidth = rhs.m_fixWidth;
m_truncatedRMS = rhs.m_truncatedRMS;
m_fractionRMS = rhs.m_fractionRMS;
@@ -205,13 +207,12 @@ InDetBeamSpotVertex::InDetBeamSpotVertex( const InDetBeamSpotVertex& rhs) :
}
StatusCode InDetBeamSpotVertex::initialize() {
- msg().setLevel( outputLevel() );
- msg(MSG::DEBUG) << "In initialize" << endreq;
+ ATH_MSG_DEBUG( "In initialize" );
return StatusCode::SUCCESS;
}
StatusCode InDetBeamSpotVertex::finalize() {
- msg(MSG::DEBUG) << "In finalize" << endreq;
+ ATH_MSG_DEBUG( "In finalize" );
clear(); // clear the data;
return StatusCode::SUCCESS;
@@ -253,8 +254,7 @@ InDetBeamSpotVertex::FitStatus InDetBeamSpotVertex::fit(std::vector< BeamSpot::V
m_getLLres = true; //set system to use these results
llResult = solveLL();
if (!llResult) {
- msg(MSG::WARNING) << "Log-likelihood fit failed: Reverting to chi2 solution"
- << endreq;
+ ATH_MSG_WARNING( "Log-likelihood fit failed: Reverting to chi2 solution" );
m_getLLres = false;
chi2Result = solveChi2(); //unfortunately need to resolve, to set correct fit parameters
}
@@ -264,8 +264,7 @@ InDetBeamSpotVertex::FitStatus InDetBeamSpotVertex::fit(std::vector< BeamSpot::V
m_getLLres = true; //set system to use these results
llResult = solveLL();
if (!llResult) {
- msg(MSG::WARNING) << "Log-likelihood fit failed: Reverting to chi2 solution"
- << endreq;
+ ATH_MSG_WARNING( "Log-likelihood fit failed: Reverting to chi2 solution" );
m_getLLres = false;
chi2Result = solveChi2(); //unfortunately need to resolve, to set correct fit parameters
}
@@ -279,32 +278,32 @@ InDetBeamSpotVertex::FitStatus InDetBeamSpotVertex::fit(std::vector< BeamSpot::V
m_fitStatus = (chi2Result ? successful : failed);
for ( std::vector< BeamSpot::VrtHolder >::iterator it =
- m_vertexData.begin();
- it != m_vertexData.end(); ++it) {
+ m_vertexData.begin();
+ it != m_vertexData.end(); ++it) {
if ( it->valid ) m_nUsed++;
}
if (msgLvl(MSG::INFO)) {
//display results
- msg(MSG::INFO) << "Fit result was: " << ( m_fitStatus == successful ? "Successful" : "Failure")
- << " using Vertex " << (m_getLLres ? "Log-likelihood":"Chi2") << " method" <<endreq;
- msg(MSG::INFO) << "Number of vertices: " << m_vertexCount << endreq;
- msg(MSG::INFO) << "x(z=0) = " << getX(0.) << " +- " << getErrX(0) << endreq;
- msg(MSG::INFO) << "y(z=0) = " << getY(0.) << " +- " << getErrY(0) << endreq;
- msg(MSG::INFO) << "z = " << getZ() << " +- " << getErrZ() << endreq;
+ ATH_MSG_INFO( "Fit result was: " << ( m_fitStatus == successful ? "Successful" : "Failure")
+ << " using Vertex " << (m_getLLres ? "Log-likelihood":"Chi2") << " method" );
+ ATH_MSG_INFO( "Number of vertices: " << m_vertexCount );
+ ATH_MSG_INFO( "x(z=0) = " << getX(0.) << " +- " << getErrX(0) );
+ ATH_MSG_INFO( "y(z=0) = " << getY(0.) << " +- " << getErrY(0) );
+ ATH_MSG_INFO( "z = " << getZ() << " +- " << getErrZ() );
double z = getZ();
- msg(MSG::INFO) << "x(z) = " << getX(z) << " +- " << getErrX(z) << endreq;
- msg(MSG::INFO) << "y(z) = " << getY(z) << " +- " << getErrY(z) << endreq;
+ ATH_MSG_INFO( "x(z) = " << getX(z) << " +- " << getErrX(z) );
+ ATH_MSG_INFO( "y(z) = " << getY(z) << " +- " << getErrY(z) );
- msg(MSG::INFO) << "TiltX = " << getTiltX() << " +- " << getErrTiltX() << endreq;
- msg(MSG::INFO) << "TiltY = " << getTiltY() << " +- " << getErrTiltY() << endreq;
+ ATH_MSG_INFO( "TiltX = " << getTiltX() << " +- " << getErrTiltX() );
+ ATH_MSG_INFO( "TiltY = " << getTiltY() << " +- " << getErrTiltY() );
- msg(MSG::INFO) << "SigmaX(z=0) = " << getSigmaX(0.) << " +- " << getErrSigmaX(0.) << endreq;
- msg(MSG::INFO) << "SigmaY(z=0) = " << getSigmaY(0.) << " +- " << getErrSigmaY(0.) << endreq;
- msg(MSG::INFO) << "SigmaZ = " << getSigmaZ() << " +- " << getErrSigmaZ() << endreq;
+ ATH_MSG_INFO( "SigmaX(z=0) = " << getSigmaX(0.) << " +- " << getErrSigmaX(0.) );
+ ATH_MSG_INFO( "SigmaY(z=0) = " << getSigmaY(0.) << " +- " << getErrSigmaY(0.) );
+ ATH_MSG_INFO( "SigmaZ = " << getSigmaZ() << " +- " << getErrSigmaZ() );
- msg(MSG::INFO) << "rhoXY = " << getRhoXY() << " +- " << getErrRhoXY() << endreq;
- msg(MSG::INFO) << "K = " << getK() << " +- " << getErrK() << endreq;
+ ATH_MSG_INFO( "rhoXY = " << getRhoXY() << " +- " << getErrRhoXY() );
+ ATH_MSG_INFO( "K = " << getK() << " +- " << getErrK() );
}
return m_fitStatus;
@@ -397,8 +396,8 @@ bool InDetBeamSpotVertex::solveLL() {
doFit2(minuit);
successfulFit(minuit, status);
if ( status.first == 3 && ( status.second == "SUCCESSFUL" ||
- status.second == "CONVERGED " ||
- status.second == "CONVERGED") ){
+ status.second == "CONVERGED " ||
+ status.second == "CONVERGED") ){
goodFit = true;
} else {
goodFit = false;
@@ -406,10 +405,10 @@ bool InDetBeamSpotVertex::solveLL() {
}
if (goodFit){
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) << "Successful fit" << endreq;
+ ATH_MSG_INFO( "Successful fit" );
setOutput( minuit);
} else {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "No LL fit convergence" << endreq;
+ ATH_MSG_DEBUG( "No LL fit convergence" );
}
delete minuit;
@@ -480,12 +479,11 @@ int InDetBeamSpotVertex::setInitialPars( TMinuit * minuit) {
}
bool InDetBeamSpotVertex::solveChi2() {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Attempting to solve " << endreq;
+ ATH_MSG_DEBUG( "Attempting to solve " );
//calulate a solution for the current set of data
if (m_vertexCount == 0) {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "No vertices stored for solve"
- << endreq;
+ ATH_MSG_DEBUG( "No vertices stored for solve" );
m_fitStatus = failed;
return false;
}
@@ -494,35 +492,33 @@ bool InDetBeamSpotVertex::solveChi2() {
//invert matrix
int fail;
m_V = m_cov;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Matrix prior to inversion\n" << m_V<< endreq;
+ ATH_MSG_DEBUG( "Matrix prior to inversion\n" << m_V);
m_V.invert(fail); //replaces the matrix with inverse
if (fail != 0) {
return false;
}
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Matrix post inversion\n" << m_V<< endreq;
+ ATH_MSG_DEBUG( "Matrix post inversion\n" << m_V);
m_p = m_V*m_x; // calculate parameter estimates
m_zSolved = m_z / m_zErr; // weighted
m_zErrSolved = 1./sqrt(m_zErr);
- if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << "Fit solved succesfully from "
- << m_vertexCount << " vertices"<< endreq;
- msg(MSG::DEBUG) << "Chi2 Solution:\n"<<m_p <<"\n"
- << "Covariance: \n" << m_V << endreq;
- }
+
+ ATH_MSG_DEBUG( "Fit solved succesfully from " << m_vertexCount << " vertices");
+ ATH_MSG_DEBUG( "Chi2 Solution:\n"<<m_p <<"\n" << "Covariance: \n" << m_V );
+
return true;
}
bool InDetBeamSpotVertex::applyOutlierRemoval() {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In Outlier removal" << endreq;
+ ATH_MSG_DEBUG( "In Outlier removal" );
// apply selection using results of a chi2 fit (eg. displacement)
if (!m_vertexData.size()) {
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) << "No vertices found" << endreq;
+ ATH_MSG_INFO( "No vertices found" );
return false;
}
static int rCount(0); // counter used for recussive mode
@@ -538,8 +534,8 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
int count(0);
for ( std::vector< BeamSpot::VrtHolder >::iterator it =
- m_vertexData.begin();
- it != m_vertexData.end(); ++it) {
+ m_vertexData.begin();
+ it != m_vertexData.end(); ++it) {
if (!it->valid) continue;
if (!m_truncatedRMS) {
@@ -568,13 +564,13 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
double mediany = (vy.size() > 1 ? vy.at(vy.size()/2) : 0.);
double medianz = (vz.size() > 1 ? vz.at(vz.size()/2) : 0.);
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Median: x: " << medianx
- << " y: " << mediany
- << " z: " << medianz << endreq;
+ ATH_MSG_DEBUG( "Median: x: " << medianx
+ << " y: " << mediany
+ << " z: " << medianz );
if (m_truncatedRMS) {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Truncating RMS ... " << endreq;
+ ATH_MSG_DEBUG( "Truncating RMS ... " );
// Only use the fraction (95%) of tracks closest to the median position to determin the RMS
int nvtx(vx.size());
std::sort(vx.begin(), vx.end(), SortDistToMedian(medianx));
@@ -586,7 +582,7 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
std::sort(vz.begin(), vz.end(), SortDistToMedian(medianz));
vz.erase(vz.begin()+int(m_fractionRMS*nvtx), vz.end());
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "... using " << int(m_fractionRMS*nvtx) << " from " << nvtx << " vertices" << endreq;
+ ATH_MSG_DEBUG( "... using " << int(m_fractionRMS*nvtx) << " from " << nvtx << " vertices" );
for (unsigned int ivtx(0); ivtx < vx.size(); ++ivtx) {
double x = vx.at(ivtx);
@@ -621,13 +617,13 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
rmsY = m_init_sy;
// should we include z here too??
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Setting RMS in x to " << rmsX << endreq;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Setting RMS in y to " << rmsY << endreq;
+ ATH_MSG_DEBUG( "Setting RMS in x to " << rmsX );
+ ATH_MSG_DEBUG( "Setting RMS in y to " << rmsY );
}
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "mean, RMS; x: " << meanx << " " << rmsX
- << " " << ", y: " << meany << " " << rmsY << ", z: " << meanz << " " << rmsZ << endreq;
+ ATH_MSG_DEBUG( "mean, RMS; x: " << meanx << " " << rmsX
+ << " " << ", y: " << meany << " " << rmsY << ", z: " << meanz << " " << rmsZ );
// varivables to hold the new chi2 and weighted z-mean values
@@ -640,8 +636,8 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
// also redetermine a new chi2 value, based on passed vertices
long failCount(0);
for ( std::vector< BeamSpot::VrtHolder >::iterator it =
- m_vertexData.begin();
- it != m_vertexData.end(); ++it) {
+ m_vertexData.begin();
+ it != m_vertexData.end(); ++it) {
if (!it->valid) continue;
int fail=0;
//if ( fabs( meanx - it->x ) > m_sigTr * rmsX ) fail += 1;
@@ -652,27 +648,27 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
/*
if ( (meanx - it->x)*(meanx-it->x)/rmsX/rmsX + (meany-it->y)*(meany-it->y)/rmsY/rmsY > m_sigTr*m_sigTr) {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Vertex info: extended past radial extent: sig."
- << sqrt((meanx - it->x)*(meanx-it->x)/rmsX/rmsX + (meany-it->y)*(meany-it->y)/rmsY/rmsY) << " > "
- << sqrt( m_sigTr*m_sigTr ) << "." << endreq;
+ if (msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG( "Vertex info: extended past radial extent: sig."
+ << sqrt((meanx - it->x)*(meanx-it->x)/rmsX/rmsX + (meany-it->y)*(meany-it->y)/rmsY/rmsY) << " > "
+ << sqrt( m_sigTr*m_sigTr ) << "." );
// don't remove these vertices yet
}
*/
if ( (medianx - it->x)*(medianx-it->x)/rmsX/rmsX + (mediany-it->y)*(mediany-it->y)/rmsY/rmsY > m_sigTr*m_sigTr) {
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Vertex info: extended past radial extent: sig."
- << sqrt((medianx - it->x)*(medianx-it->x)/rmsX/rmsX + (mediany-it->y)*(mediany-it->y)/rmsY/rmsY) << " > "
- << sqrt( m_sigTr*m_sigTr ) << "." << endreq;
+ ATH_MSG_DEBUG( "Vertex info: extended past radial extent: sig."
+ << sqrt((medianx - it->x)*(medianx-it->x)/rmsX/rmsX + (mediany-it->y)*(mediany-it->y)/rmsY/rmsY) << " > "
+ << sqrt( m_sigTr*m_sigTr ) << "." );
fail += 128;
}
//and in Z?
if (fail) { // TBD only allow failed vertices here to be removed on every nth iteration (to aid stability)
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Vertex reject from simple mean; reason: " << fail << " : x,y,z: "
+ ATH_MSG_DEBUG( "Vertex reject from simple mean; reason: " << fail << " : x,y,z: "
<< it->x << " " << it->y << " " << it->z
<< " , sigma(x,y,z): " << sqrt(it->vxx) << " " << sqrt(it->vyy)
<< " " << sqrt(it->vzz)
- << endreq;
+ );
it->valid = false;
// it->outlierRemoved = true;
@@ -702,7 +698,7 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
}// end of valid vertex
}// for
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Removed: " << failCount << " vertices from simple mean,RMS." << endreq;
+ ATH_MSG_DEBUG( "Removed: " << failCount << " vertices from simple mean,RMS." );
// update the new global chi2 - need to do this after the invert has been successful though !!!
@@ -722,12 +718,12 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
// make comparisons of old chi2, new mean and new chi2 values
if (msgLvl(MSG::DEBUG) ) {
- msg(MSG::DEBUG) << "Mean position: x,y,z " << meanx << " " << meany << " " << meanz << endreq;
- msg(MSG::DEBUG) << " RMS: x,y,z " << rmsX << " " << rmsY << " " << rmsZ << endreq;
+ ATH_MSG_DEBUG( "Mean position: x,y,z " << meanx << " " << meany << " " << meanz );
+ ATH_MSG_DEBUG( " RMS: x,y,z " << rmsX << " " << rmsY << " " << rmsZ );
- msg(MSG::DEBUG) << "Original chi2:" << m_p << "\n" << m_V << "\n" << m_zSolved << " " << m_zErrSolved << endreq;
+ ATH_MSG_DEBUG( "Original chi2:" << m_p << "\n" << m_V << "\n" << m_zSolved << " " << m_zErrSolved );
- msg(MSG::DEBUG) << "New chi2:" << chi2Pos << "\n" << chi2Cov << "\n" << zpos << " " << zerr << endreq;
+ ATH_MSG_DEBUG( "New chi2:" << chi2Pos << "\n" << chi2Cov << "\n" << zpos << " " << zerr );
} // debug statement
@@ -739,15 +735,15 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
/* TBD
if (goodChi2) {
- if (msgLvl(MSG::DEBUG)) ) msg(MSG::DEBUG) << "Replacing original chi2 with new value" << endreq;
+ if (msgLvl(MSG::DEBUG)) ) ATH_MSG_DEBUG( "Replacing original chi2 with new value" );
}
*/
//Dubious - !!! TBD: move this into the general LL solve area and do it properly!
- m_init_sx = rmsX;
- m_init_sy = rmsY;
+ //m_init_sx = rmsX;
+ //m_init_sy = rmsY;
m_init_sz = rmsZ;
// perform LL fit?
@@ -759,7 +755,7 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
}
if ( llSolve and rCount > 0 ) {
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) << "Log-Likelihood fit converged in outlier removal. Exiting outlier removal." << endreq;
+ ATH_MSG_INFO( "Log-Likelihood fit converged in outlier removal. Exiting outlier removal." );
return true;
}
@@ -775,9 +771,9 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
// ll solution not used, or not trusted
// set a wide solution
m_getLLres = false;
- if (msgLvl(MSG::INFO)) msg(MSG::INFO) << ": removeOutliers: LL fit not use/converged/trusted - " <<
- "using chi2 for mean and simple RMS for width values " << endreq;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << llSolve << " " << getSigmaX(0.) << " " << getSigmaY(0.) << " " << getRhoXY() << endreq;
+ ATH_MSG_INFO( ": removeOutliers: LL fit not use/converged/trusted - " <<
+ "using chi2 for mean and simple RMS for width values " );
+ ATH_MSG_DEBUG( llSolve << " " << getSigmaX(0.) << " " << getSigmaY(0.) << " " << getRhoXY() );
xbar = getX(0.);
ybar = getY(0.);
ax = getTiltX();
@@ -858,24 +854,22 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
if ( !vit->second) continue; // no valid pointer ...
if ( !vit->second->valid) continue; // ignore vertices set to invalid - shouldn't happen here though
if ( fCount >= m_singleIterationMax && m_singleIterationMax != -1) {
- if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << " removeOutlier: Reached max number of vertex rejections for this iteration.\n"
- << "\tNeed to recalculate mean positions." << endreq;
- }
+ ATH_MSG_DEBUG( " removeOutlier: Reached max number of vertex rejections for this iteration.\n"
+ << "\tNeed to recalculate mean positions." );
break;
} // if reached max iterations for this round
- if (msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Vertex chi2: " << vit->first << endreq;
+ ATH_MSG_VERBOSE( "Vertex chi2: " << vit->first );
if ( vit->first < m_outlierChi2Tr ) {
if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << " No more 'bad' vertices found in this iteration." ;
+ ATH_MSG_DEBUG( " No more 'bad' vertices found in this iteration." );
if (fCount == 0) {
- msg(MSG::DEBUG) << " No futher vertices removed - moving to final iteration" << endreq;
+ ATH_MSG_DEBUG( " No futher vertices removed - moving to final iteration" );
} else {
- msg(MSG::DEBUG) << " Moving to next iteration of outlier removal." << endreq;
+ ATH_MSG_DEBUG( " Moving to next iteration of outlier removal." );
}
}
break;
@@ -885,20 +879,18 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
++fCount;
vit->second->valid = false;
// vit->second->outlierRemoved = true;
- if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << "Vertex rejected; chi2: " << vit->first <<". pos(x,y,z): "
+ ATH_MSG_DEBUG( "Vertex rejected; chi2: " << vit->first <<". pos(x,y,z): "
<< vit->second->x << " " << vit->second->y << " " << vit->second->z
<< " , sigma(x,y,z): " << sqrt(vit->second->vxx) << " " << sqrt(vit->second->vyy)
<< " " << sqrt(vit->second->vzz)
- << endreq;
- } // debug
+ );
} // for
// if no vertices removed and ll fit still fails, then we continue to have a problem ...
if (fCount == 0 && m_useLL && !llSolve && rCount !=0 ) { // if first iteration, we have another iteration later.
- // if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "In LL test MODE !!!!" << endreq;
- if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "No vertices removed and fit still fails - most likely final result will fail" << endreq;
+ // if (msgLvl(MSG::WARNING)) ATH_MSG_WARNING( "In LL test MODE !!!!" );
+ ATH_MSG_WARNING( "No vertices removed and fit still fails - most likely final result will fail" );
// this is our 'last-ditch approach'. Split the collection of vertices into two 'random' sets and solve for each.
// if both successful, then compare (and be a little confused). If only one, compare to chi2 and take if ok.
@@ -918,61 +910,61 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
bool goodFit1(false), goodFit2(false);
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Attempting fit with vextex half: 1"<< endreq;
+ ATH_MSG_DEBUG( "Attempting fit with vextex half: 1");
m_vertexData = vertex1;
llSolve = solveLL();
//m_getLLres = llSolve; // allow the log-likelihood accessor values to returned, if sucessful
m_getLLres = true;
- if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Fit using \"vertex1\" " << ( llSolve ? "Successful": "Failed") << endreq;
+ ATH_MSG_WARNING( "Fit using \"vertex1\" " << ( llSolve ? "Successful": "Failed") );
if (llSolve) {
goodFit1 = true;
if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << "x: " << getX(0) << " +- " << getErrX(0) << endreq;
- msg(MSG::DEBUG) << "y: " << getY(0) << " +- " << getErrY(0) << endreq;
- msg(MSG::DEBUG) << "z: " << getZ() << " +- " << getErrZ() << endreq;
- msg(MSG::DEBUG) << "sx: " << getSigmaX(0) << " +- " << getErrSigmaX(0) << endreq;
- msg(MSG::DEBUG) << "sy: " << getSigmaY(0) << " +- " << getErrSigmaY(0) << endreq;
- msg(MSG::DEBUG) << "sz: " << getSigmaZ() << " +- " << getErrSigmaZ() << endreq;
- msg(MSG::DEBUG) << "ax: " << getTiltX() << " +- " << getErrTiltX() << endreq;
- msg(MSG::DEBUG) << "ay: " << getTiltY() << " +- " << getErrTiltY() << endreq;
- msg(MSG::DEBUG) << "sxy:" << getSigmaXY(0) << " +- " << getErrSigmaXY(0) << endreq;
- msg(MSG::DEBUG) << "rh: " << getRhoXY() << " +- " << getErrRhoXY() << endreq;
- msg(MSG::DEBUG) << "k: " << getK() << " +- " << getErrK() << endreq;
+ ATH_MSG_DEBUG( "x: " << getX(0) << " +- " << getErrX(0) );
+ ATH_MSG_DEBUG( "y: " << getY(0) << " +- " << getErrY(0) );
+ ATH_MSG_DEBUG( "z: " << getZ() << " +- " << getErrZ() );
+ ATH_MSG_DEBUG( "sx: " << getSigmaX(0) << " +- " << getErrSigmaX(0) );
+ ATH_MSG_DEBUG( "sy: " << getSigmaY(0) << " +- " << getErrSigmaY(0) );
+ ATH_MSG_DEBUG( "sz: " << getSigmaZ() << " +- " << getErrSigmaZ() );
+ ATH_MSG_DEBUG( "ax: " << getTiltX() << " +- " << getErrTiltX() );
+ ATH_MSG_DEBUG( "ay: " << getTiltY() << " +- " << getErrTiltY() );
+ ATH_MSG_DEBUG( "sxy:" << getSigmaXY(0) << " +- " << getErrSigmaXY(0) );
+ ATH_MSG_DEBUG( "rh: " << getRhoXY() << " +- " << getErrRhoXY() );
+ ATH_MSG_DEBUG( "k: " << getK() << " +- " << getErrK() );
}// debug
} // good fit
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Attempting fit with vextex half: 2"<< endreq;
+ ATH_MSG_DEBUG( "Attempting fit with vextex half: 2");
m_vertexData = vertex2;
llSolve = solveLL();
// m_getLLres = llSolve; // allow the log-likelihood accessor values to returned, if sucessful
m_getLLres = true;
- if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Fit using \"vertex2\" " << ( llSolve ? "Successful": "Failed") << endreq;
+ ATH_MSG_WARNING( "Fit using \"vertex2\" " << ( llSolve ? "Successful": "Failed") );
if (llSolve) {
goodFit2 = true;
if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << "x: " << getX(0) << " +- " << getErrX(0) << endreq;
- msg(MSG::DEBUG) << "y: " << getY(0) << " +- " << getErrY(0) << endreq;
- msg(MSG::DEBUG) << "z: " << getZ() << " +- " << getErrZ() << endreq;
- msg(MSG::DEBUG) << "sx: " << getSigmaX(0) << " +- " << getErrSigmaX(0) << endreq;
- msg(MSG::DEBUG) << "sy: " << getSigmaY(0) << " +- " << getErrSigmaY(0) << endreq;
- msg(MSG::DEBUG) << "sz: " << getSigmaZ() << " +- " << getErrSigmaZ() << endreq;
- msg(MSG::DEBUG) << "ax: " << getTiltX() << " +- " << getErrTiltX() << endreq;
- msg(MSG::DEBUG) << "ay: " << getTiltY() << " +- " << getErrTiltY() << endreq;
- msg(MSG::DEBUG) << "sxy:" << getSigmaXY(0) << " +- " << getErrSigmaXY(0) << endreq;
- msg(MSG::DEBUG) << "rh: " << getRhoXY() << " +- " << getErrRhoXY() << endreq;
- msg(MSG::DEBUG) << "k: " << getK() << " +- " << getErrK() << endreq;
+ ATH_MSG_DEBUG( "x: " << getX(0) << " +- " << getErrX(0) );
+ ATH_MSG_DEBUG( "y: " << getY(0) << " +- " << getErrY(0) );
+ ATH_MSG_DEBUG( "z: " << getZ() << " +- " << getErrZ() );
+ ATH_MSG_DEBUG( "sx: " << getSigmaX(0) << " +- " << getErrSigmaX(0) );
+ ATH_MSG_DEBUG( "sy: " << getSigmaY(0) << " +- " << getErrSigmaY(0) );
+ ATH_MSG_DEBUG( "sz: " << getSigmaZ() << " +- " << getErrSigmaZ() );
+ ATH_MSG_DEBUG( "ax: " << getTiltX() << " +- " << getErrTiltX() );
+ ATH_MSG_DEBUG( "ay: " << getTiltY() << " +- " << getErrTiltY() );
+ ATH_MSG_DEBUG( "sxy:" << getSigmaXY(0) << " +- " << getErrSigmaXY(0) );
+ ATH_MSG_DEBUG( "rh: " << getRhoXY() << " +- " << getErrRhoXY() );
+ ATH_MSG_DEBUG( "k: " << getK() << " +- " << getErrK() );
}// debug
}
// what now ? ...
- if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Fit was " << ( goodFit2 || goodFit1 ? "Successful ": "Unsuccessful ")
- <<" using a subset of the available vertices" << endreq;
- if (msgLvl(MSG::WARNING) && ( goodFit2 || goodFit1) )
- msg(MSG::WARNING) << "Using these subset data for final result!!! " << endreq;
- if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "FIT HALFVERTX" << endreq;
+ ATH_MSG_WARNING( "Fit was " << ( goodFit2 || goodFit1 ? "Successful ": "Unsuccessful ")
+ <<" using a subset of the available vertices" );
+ if (( goodFit2 || goodFit1) )
+ ATH_MSG_WARNING( "Using these subset data for final result!!! " );
+ ATH_MSG_WARNING( "FIT HALFVERTX" );
if (goodFit1) {
m_vertexData = vertex1;
@@ -990,20 +982,15 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
// recursive mode
- if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << " Recursive debug: Loop: " << rCount << ". Number of failed vertices: " << fCount << endreq ;
- }
+ ATH_MSG_DEBUG( " Recursive debug: Loop: " << rCount << ". Number of failed vertices: " << fCount );
+
++rCount;
if ( fCount > 0 || ( fCount == 0 && rCount == 1 && !llSolve)) { // if failed vertices or, no failed, first iteration, and no succesful fit
if ( rCount > m_maxOutlierLoops) {
- if (msgLvl(MSG::WARNING)) {
- msg(MSG::WARNING) << "OutlierRemoval: Reached maximum number of recursive loops: " << rCount
- << ". No more iterations performed." << endreq;
- }
+ ATH_MSG_WARNING( "OutlierRemoval: Reached maximum number of recursive loops: " << rCount
+ << ". No more iterations performed." );
} else {
- if (msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG) << "OutlierRemoval: Entering recursive loop: " << rCount << endreq;
- }
+ ATH_MSG_DEBUG( "OutlierRemoval: Entering recursive loop: " << rCount );
applyOutlierRemoval();
} // if entering loop
} // if fails > 0
@@ -1013,7 +1000,7 @@ bool InDetBeamSpotVertex::applyOutlierRemoval() {
bool InDetBeamSpotVertex::successfulFit( TMinuit * minuit,
- std::pair<int, std::string> & status) {
+ std::pair<int, std::string> & status) {
if (!minuit) return false;
//This should be called directly after the fit
std::string sRes = minuit->fCstatu.Data();
@@ -1022,8 +1009,7 @@ bool InDetBeamSpotVertex::successfulFit( TMinuit * minuit,
Int_t npari,nparx,istat;
minuit->mnstat(fmin, fedm, errdef,npari,nparx,istat);
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit reports status: " << istat << " and "
- << sRes << endreq;
+ ATH_MSG_DEBUG( "Fit reports status: " << istat << " and " << sRes );
status.first = istat;
status.second = sRes;
@@ -1034,25 +1020,24 @@ bool InDetBeamSpotVertex::successfulFit( TMinuit * minuit,
minuit->GetParameter(6,x,ex); // rhoxy
if ( fabs(x) > m_rhoFail){
sanityPassed = false;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed with rhoxy: " << x << " > " << m_rhoFail << endreq;
+ ATH_MSG_DEBUG( "Fit Failed with rhoxy: " << x << " > " << m_rhoFail );
}
minuit->GetParameter(4,x,ex); // sigma x
if ( fabs(x) < m_widthFail ){
sanityPassed = false;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed with sigmaX:" << x << " > " << m_widthFail << endreq;
-
+ ATH_MSG_DEBUG( "Fit Failed with sigmaX:" << x << " > " << m_widthFail );
}
minuit->GetParameter(5,x,ex); // sigma y
if ( fabs(x) < m_widthFail ){
sanityPassed = false;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed with sigmaY: " << x << " > " <<m_widthFail << endreq;
+ ATH_MSG_DEBUG( "Fit Failed with sigmaY: " << x << " > " <<m_widthFail );
}
minuit->GetParameter(7,x,ex); // k
if ( fabs(x) < m_kMinFail || fabs(x) > m_kMaxFail ){
sanityPassed = false;
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed with k: " << x << " > " << m_kMaxFail
- << ", or " << x << " < " << m_kMinFail << endreq;
+ ATH_MSG_DEBUG( "Fit Failed with k: " << x << " > " << m_kMaxFail
+ << ", or " << x << " < " << m_kMinFail );
}
@@ -1062,7 +1047,7 @@ bool InDetBeamSpotVertex::successfulFit( TMinuit * minuit,
status.first = 99;
status.second = "FAILED BEAMSPOT SANITY CHECK";
}
- if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit " << ( sanityPassed ? "Passed": "Failed") << " sanity check: " << endreq;
+ ATH_MSG_DEBUG( "Fit " << ( sanityPassed ? "Passed": "Failed") << " sanity check: " );
if ( istat != 3) return false;
@@ -1148,10 +1133,10 @@ double BeamSpot::pdfxy(double *, double *) {
double pxy = 1/(2*Pi*det)*exp( -0.5* ( ( x - mux -ax*z) * covXX * ( x - mux - ax*z)
- + ( y - muy - ay*z) * covYY * ( y - muy - ay*z)
- + 2*( x - mux - *z) * covXY * ( y - muy - ay*z)
- )
- );
+ + ( y - muy - ay*z) * covYY * ( y - muy - ay*z)
+ + 2*( x - mux - *z) * covXY * ( y - muy - ay*z)
+ )
+ );
return pxy;
*/
return 0; // TBD dlete
@@ -1222,8 +1207,8 @@ void BeamSpot::myFCN_LLsolverNorm( Int_t &, Double_t *, Double_t & /*f*/, Double
+ ( y - par[1] - par[3]*z) * covYY * ( y - par[1] - par[3]*z)
+ 2*( x - par[0] - par[2]*z) * covXY * ( y - par[1] - par[3]*z)
);
-
- temp += TMath::Log( 2*Pi * par[9]*par[9] ) + ( z - par[8]) * (z-par[8]) / (par[9] * par[9] );
+
+ temp += TMath::Log( 2*Pi * par[9]*par[9] ) + ( z - par[8]) * (z-par[8]) / (par[9] * par[9] );
*/
// apply normalisations
@@ -1257,6 +1242,11 @@ void InDetBeamSpotVertex::doFit2( TMinuit * minuit) {
minuit->FixParameter(2);
minuit->FixParameter(3);
+ if(m_fixWidth){
+ minuit->FixParameter(4);
+ minuit->FixParameter(5);
+ }
+
minuit->Migrad();
minuit->Migrad();
//release k,rho
@@ -1265,9 +1255,10 @@ void InDetBeamSpotVertex::doFit2( TMinuit * minuit) {
minuit->Migrad();
}
- minuit->Release(6);
- minuit->Migrad();
-
+ if(!m_fixWidth){
+ minuit->Release(6);
+ minuit->Migrad();
+ }
minuit->Release(0);
minuit->Release(1);
@@ -1293,37 +1284,63 @@ std::map<std::string,double> InDetBeamSpotVertex::getCovMap() const {
//This is the method that was used before to put the covariance matrix in the required order
//We don't want to mess with this, because no one knows the original order
- static const int map[] = {0,0,6,-1,-1,-1,-1,2,-2,-2};
+ //static const int map[] = {0,0,6,-1,-1,-1,-1,2,-2,-2};
+ int map[] = {1,2,9,3,4,5,6,10,7,8};
+ if(m_fixInputK){
+ // 1 2 3 4 5 6 7 8 9 10
+ //int map2[] = {0,0,5,-1,-1,-1,-1,0,-2,-2};
+ int map2[] = {1,2,8,3,4,5,6,9,7,10};
+ for(int i=0; i < 10; ++i){
+ map[i] = map2[i];
+ }
+ } else if (m_fixWidth) {
+ //int map2[] = {0,0,3,-1,-1,-5,-6,-1, -8,-5};
+ int map2[] = {1,2,6,3,4,8,9,7,10,5};
+ for(int i=0; i < 10; ++i){
+ map[i] = map2[i];
+ }
+ }
+
int temp = 0;
for (int i=0;i<10;++i) {
for (int j=i;j<10;++j) {
- covVector[temp++] = m_VLL(i + 1 +map[i], j+1+ map[j]);
+ if( m_fixInputK && (i == 9 || j ==9 )){
+ covVector[temp++] = 0;
+ } else if ( m_fixWidth && ( i == 5 || i == 6 || i == 8 || j == 5 || j == 6 || j == 8 ) ){
+ covVector[temp++] = 0;
+ }else{
+ covVector[temp++] = m_VLL( map[i], map[j] );
+ }
}
}
//This array is in the order required from the original ntuple format
- const std::string keyArr[] = {"posXErr","covXY","covXZ","covXTiltX","covXTiltY","covXSx","covXSy","covXSz",
- "covXRhoXY","covXk","posYErr","covYZ","covYTiltX","covYTiltY","covYSx","covYSy",
- "covYSz","covYRhoXY","covYk","posZErr","covZTiltX","covZTiltY","covZSx","covZSy","covZSz",
- "covZRhoXY","covZk","tiltXErr","covTiltXTiltY","covTiltXSx","covTiltXSy","covTiltXSz","covTiltXRhoXY",
- "covTiltXk","tiltYErr","covTiltYSx","covTiltYSy","covTiltYSz","covTiltYRhoXY","covTiltYk","sigmaXErr",
- "covSxSy","covSxSz","covSxRhoXY","covSxk","sigmaYErr","covSySz","covSyRhoXY","covSyk","sigmaZErr","covSzRhoXY",
- "covSzk","rhoXYErr","covRhoXYk","kErr"};
+ const std::string keyArr[] = {"posXErr","covXY","covXZ","covXTiltX","covXTiltY","covXSx","covXSy","covXSz","covXRhoXY","covXk",
+ "posYErr","covYZ","covYTiltX","covYTiltY","covYSx","covYSy","covYSz","covYRhoXY","covYk",
+ "posZErr","covZTiltX","covZTiltY","covZSx","covZSy","covZSz","covZRhoXY","covZk",
+ "tiltXErr","covTiltXTiltY","covTiltXSx","covTiltXSy","covTiltXSz","covTiltXRhoXY","covTiltXk",
+ "tiltYErr","covTiltYSx","covTiltYSy","covTiltYSz","covTiltYRhoXY","covTiltYk",
+ "sigmaXErr","covSxSy","covSxSz","covSxRhoXY","covSxk",
+ "sigmaYErr","covSySz","covSyRhoXY","covSyk",
+ "sigmaZErr","covSzRhoXY","covSzk",
+ "rhoXYErr","covRhoXYk",
+ "kErr"};
/* const std::string keyArr2[] = {"posXErr","covYX","covZX","covTiltXX","covTiltYX","covSxX","covSyX","covSzX",
- "covRhoXYX","covkX","posYErr","covZY","covTiltXY","covTiltYY","covSxY","covSyY",
- "covSzY","covRhoXYY","covkY","posZErr","covTiltXZ","covTiltYZ","covSxZ","covSyZ","covSzZ",
- "covRhoXYZ","covkZ","tiltXErr","covTiltYTiltX","covSxTiltX","covSyTiltX","covSzTiltX","covRhoXYTiltX",
- "covkTiltX","tiltYErr","covSxTiltY","covSyTiltY","covSzTiltY","covRhoXYTiltY","covkTiltY","sigmaXErr",
- "covSySx","covSzSx","covRhoXYSx","covkSx","sigmaYErr","covSzSy","covRhoXYSy","covkSy","sigmaZErr","covRhoXYSz",
- "covkSz","rhoXYErr","covkRhoXY","kErr"};
+ "covRhoXYX","covkX","posYErr","covZY","covTiltXY","covTiltYY","covSxY","covSyY",
+ "covSzY","covRhoXYY","covkY","posZErr","covTiltXZ","covTiltYZ","covSxZ","covSyZ","covSzZ",
+ "covRhoXYZ","covkZ","tiltXErr","covTiltYTiltX","covSxTiltX","covSyTiltX","covSzTiltX","covRhoXYTiltX",
+ "covkTiltX","tiltYErr","covSxTiltY","covSyTiltY","covSzTiltY","covRhoXYTiltY","covkTiltY","sigmaXErr",
+ "covSySx","covSzSx","covRhoXYSx","covkSx","sigmaYErr","covSzSy","covRhoXYSy","covkSy","sigmaZErr","covRhoXYSz",
+ "covkSz","rhoXYErr","covkRhoXY","kErr"};
*/
//Now that everything should be in the right order, it's simple to set the covariance matrix map correctly:
for(int i = 0; i < 55; i++){
covMap[keyArr[i]] = covVector[i];
+ //std::cout << keyArr[i] << " " << covVector[i] << std::endl;
//covMap[keyArr2[i]]= covVector[i];
}
@@ -1345,10 +1362,10 @@ std::map<std::string,double> InDetBeamSpotVertex::getCovMap() const {
double z = getZ();
CLHEP::HepSymMatrix covc = getCov(z);
- covMap["posXErr"] = sqrt( covc(1,1) ); //xcxc
- covMap["posYErr"] = sqrt( covc(2,2) ); //ycyc
- covMap["tiltXErr"] = sqrt( covc(3,3) ); //axcaxc
- covMap["tiltYErr"] = sqrt( covc(4,4) ); //aycayc
+ covMap["posXErr"] = sqrt( covc(1,1) ); //xcxc
+ covMap["posYErr"] = sqrt( covc(2,2) ); //ycyc
+ covMap["tiltXErr"] = sqrt( covc(3,3) ); //axcaxc
+ covMap["tiltYErr"] = sqrt( covc(4,4) ); //aycayc
covMap["sigmaXErr"] = sqrt( getErrSigmaX(z)*getErrSigmaX(z) ); //sxcsxc
covMap["sigmaYErr"] = sqrt( getErrSigmaY(z)*getErrSigmaY(z) ); //sycsyc
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.h b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.h
index 25a03777f6ab..19625b88916c 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.h
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/InDetBeamSpotVertex.h
@@ -60,7 +60,7 @@ namespace InDet {
sqrt(m_V(3,3) + 2*z*m_V(4,3) + z*z*m_V(4,4)) );
}
double getErrZ() const {
- return (m_getLLres ? sqrt(m_VLL(9,9)):m_zErrSolved);} // defined as centroid of beam
+ return (m_getLLres ? (m_fixInputK ? sqrt(m_VLL(8,8)) : sqrt(m_VLL(9,9))):m_zErrSolved);} // defined as centroid of beam
double getSigmaX(double) const {
return (m_getLLres ? m_pLL(5) : m_def_sx);}
@@ -73,7 +73,7 @@ namespace InDet {
double getErrSigmaY(double) const {
return (m_getLLres ? sqrt(m_VLL(6,6)) :0.);}
double getErrSigmaZ() const {
- return (m_getLLres ? sqrt(m_VLL(10,10)) :0.);}
+ return (m_getLLres ? (m_fixInputK ? sqrt(m_VLL(9,9)) : sqrt(m_VLL(10,10))) :0.);}
double getTiltX() const {
return (m_getLLres ? m_pLL(3) :m_p(2));}
@@ -88,7 +88,7 @@ namespace InDet {
double getRhoXY() const { return (m_getLLres ? m_pLL(7) : 0.);}
double getK() const { return (m_getLLres ? m_pLL(8) : 0.);}
double getErrRhoXY() const { return (m_getLLres ? sqrt(m_VLL(7,7)) : 0.);}
- double getErrK() const { return (m_getLLres ? sqrt(m_VLL(8,8)) : 0.);}
+ double getErrK() const { return (m_getLLres ? (m_fixInputK ? 0. : sqrt(m_VLL(8,8)) ) : 0.);}
double getSigmaXY(double z) const {return getRhoXY()*getSigmaX(z)*getSigmaY(z);}
@@ -192,6 +192,7 @@ namespace InDet {
bool m_fixInputK; // fix in LL fit the K factor to it's initial value
bool m_useLLNorm; // use the normalisation mode of the pdf...
+ bool m_fixWidth;
int m_nUsed;
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.cxx
new file mode 100644
index 000000000000..088b0389fc9f
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.cxx
@@ -0,0 +1,382 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// @file RefitTracksAndVertex.cxx
+// Anthony Morley
+//
+
+#include "RefitTracksAndVertex.h"
+#include "GaudiKernel/ListItem.h"
+#include "GaudiKernel/SystemOfUnits.h"
+#include "TrkParameters/TrackParameters.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrkRIO_OnTrack/RIO_OnTrack.h"
+#include "TrkFitterInterfaces/ITrackFitter.h"
+#include "TrkVertexFitterInterfaces/IVertexFitter.h"
+#include "TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h"
+#include "Identifier/Identifier.h"
+#include "AtlasDetDescr/AtlasDetectorID.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/VertexAuxContainer.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "VxVertex/VxCandidate.h"
+
+
+
+RefitTracksAndVertex::RefitTracksAndVertex( const std::string& name, ISvcLocator* pSvcLocator )
+ : AthAlgorithm (name, pSvcLocator)
+ , m_idHelper(nullptr)
+ , m_vertexListInput("PrimaryVertices")
+ , m_trackListOutput("SiTracksWithRefitTracksAndVertex")
+ , m_outputVertexContainerName("RefittedPrimaryVertex")
+ , m_trackFitter("Trk::KalmanFitter/TrkKalmanFitter")
+ , m_vertexFitter("Trk::FullVertexFitter")
+ , m_applyTrkSel(false)
+ , m_refitTracks(true)
+ , m_addPM(false)
+ , m_selPtMin(0)
+ , m_selEtaMin(-2.5)
+ , m_selEtaMax(2.5)
+ , m_selNHitPIXMin(0)
+ , m_selNHitSCTMin(6)
+ , m_nTracksProcessed(0)
+ , m_nTracksPresel(0)
+ , m_nTracksAccepted(0)
+ , m_nRejectPIX(0)
+ , m_nRejectSCT(0)
+ , m_nRejectTRT(0)
+ , m_nRejectPM(0)
+{
+ // declare algorithm parameters
+ declareProperty("VertexListInput" , m_vertexListInput) ;
+ declareProperty("TrackListOutput" , m_trackListOutput) ;
+ declareProperty("OutputVertexContainer" , m_outputVertexContainerName) ;
+ declareProperty("TrackFitter" , m_trackFitter) ;
+ declareProperty("VertexFitterTool" , m_vertexFitter);
+ declareProperty("ApplyTrkSel" , m_applyTrkSel) ;
+ declareProperty("RefitTracks" , m_refitTracks) ;
+ declareProperty("AddPseudoMeasurement", m_addPM);
+ declareProperty("SelPtMin" , m_selPtMin) ;
+ declareProperty("SelEtaMin" , m_selEtaMin) ;
+ declareProperty("SelEtaMax" , m_selEtaMax) ;
+ declareProperty("SelNHitPIXMin" , m_selNHitPIXMin) ;
+ declareProperty("SelNHitSCTMin" , m_selNHitSCTMin) ;
+
+}
+
+RefitTracksAndVertex::~RefitTracksAndVertex() {}
+
+StatusCode RefitTracksAndVertex::initialize() {
+ //retrieve the DetectorStore service
+ StatusCode status=detStore().retrieve() ;
+ if( status.isFailure() ) {
+ ATH_MSG_FATAL( "DetectorStore service not found !" );
+ return status;
+ }else{
+ ATH_MSG_DEBUG( "DetectorStore retrieved!" );
+ }
+
+ if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) {
+ ATH_MSG_FATAL( "Could not get AtlasDetectorID helper" );
+ return StatusCode::FAILURE;
+ }
+
+ // Set up fitter
+ status = m_trackFitter.retrieve();
+ if (status.isFailure()) {
+ ATH_MSG_FATAL( "Could not find tool " << m_trackFitter << ". Exiting.");
+ return status ;
+ } else {
+ ATH_MSG_DEBUG( " Got " << m_trackFitter << " as TrackFitter. " ) ;
+ }
+
+ // Set up fitter
+ status = m_vertexFitter.retrieve();
+ if (status.isFailure()) {
+ ATH_MSG_FATAL( "Could not find tool " << m_vertexFitter << ". Exiting.");
+ return status ;
+ } else {
+ ATH_MSG_DEBUG( " Got " << m_vertexFitter << " as VertexFitter. " ) ;
+ }
+
+ // Print input properties
+ if( m_applyTrkSel ) {
+ ATH_MSG_DEBUG(
+ "Track selection will be applied:"
+ << "\n " << m_selEtaMin << " < eta < " << m_selEtaMax
+ << "\n pT > " << m_selPtMin/Gaudi::Units::GeV << " GeV"
+ << "\n number_of_PIXEL_hits >= " << m_selNHitPIXMin
+ << "\n number_of_SCT_hits >= " << m_selNHitSCTMin
+ ) ;
+ } else {
+ ATH_MSG_DEBUG( "NO selection will be applied to tracks" ) ;
+ }
+ return StatusCode::SUCCESS;
+} // initialize(.)
+
+StatusCode RefitTracksAndVertex::finalize() {
+ ATH_MSG_DEBUG( (*this) ) ;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode RefitTracksAndVertex::execute() {
+
+ TrackCollection* outputtracks = new TrackCollection() ;
+
+ //---------------------------
+ // Primary Vertex
+ //---------------------------
+ const xAOD::VertexContainer* vertices = 0;
+ if ( !evtStore()->retrieve( vertices, m_vertexListInput ).isSuccess() ){ // retrieve arguments: container type, container key
+ ATH_MSG_WARNING("execute() Failed to retrieve Reconstructed vertex container. " << m_vertexListInput );
+ //ATH_MSG_ERROR(evtStore()->dump());
+ return StatusCode::SUCCESS;
+ }
+
+ if( evtStore()->record( outputtracks, m_trackListOutput ).isFailure() ) {
+ ATH_MSG_ERROR( "Failed to record trackcollection with name " << m_trackListOutput );
+ }
+
+ xAOD::VertexContainer* theVertexContainer = new xAOD::VertexContainer;
+ xAOD::VertexAuxContainer* theVertexAuxContainer = new xAOD::VertexAuxContainer;
+ theVertexContainer->setStore( theVertexAuxContainer );
+
+ CHECK( evtStore()->record(theVertexContainer, m_outputVertexContainerName) );
+ CHECK( evtStore()->record(theVertexAuxContainer, m_outputVertexContainerName + "Aux.") );
+
+ const xAOD::Vertex* primaryVertex = 0;
+ for( auto vertex: *vertices ) {
+ if( vertex->vertexType() == xAOD::VxType::PriVtx ) {
+ primaryVertex = vertex;
+ break;
+ }
+ }
+
+ std::vector< const Trk::TrackParameters* > trackParametersToFit;
+
+ // Check we have found a PV
+ if( primaryVertex ){
+ if(primaryVertex->nTrackParticles() < 2)
+ return StatusCode::SUCCESS;
+
+ //Loop over all tracks in the PV
+ for(unsigned int i =0; i< primaryVertex->nTrackParticles(); ++i ){
+ auto trackParticle = primaryVertex->trackParticle( i );
+ auto track = trackParticle->track();
+ // Check that there is Trk::Track
+ if(track){
+ //Refit Track using on the SCT hits
+ if(m_refitTracks){
+ auto newTrack = fitSCTOnlyTrack( track ) ;
+ if(newTrack){
+ outputtracks->push_back( newTrack ) ;
+ if(newTrack->perigeeParameters())
+ trackParametersToFit.push_back( newTrack->perigeeParameters() );
+ }
+ } else {
+ if(track->perigeeParameters())
+ trackParametersToFit.push_back( track->perigeeParameters() );
+ }
+ }
+ }
+ m_nTracksProcessed += primaryVertex->nTrackParticles() ;
+ m_nTracksAccepted += outputtracks->size() ;
+ // Make sure there are at least 2 tracks to fit in the vertex fit
+
+ ATH_MSG_DEBUG(" From: " << primaryVertex->nTrackParticles() << " tracks " << trackParametersToFit.size() << " will be used" );
+
+ if(trackParametersToFit.size() > 1){
+
+ Amg::Vector3D position = primaryVertex->position();
+ // Fit and store vertex
+ auto vertex = m_vertexFitter->fit(trackParametersToFit, position );
+ if(vertex){
+ theVertexContainer->push_back (vertex);
+ vertex->setVertexType(xAOD::VxType::PriVtx);
+
+ ATH_MSG_DEBUG("Old Vtx " << primaryVertex->position());
+ ATH_MSG_DEBUG("New Vtx " << vertex->position());
+ }
+ }
+ }
+
+
+
+
+
+ return StatusCode::SUCCESS;
+} // execute(.)
+
+
+const Trk::PseudoMeasurementOnTrack* RefitTracksAndVertex::createPMfromSi ( const Trk::Perigee* mp ) {
+ Trk::DefinedParameter z0 ( mp->parameters()[Trk::z0], Trk::z0 ) ;
+ Trk::DefinedParameter theta( mp->parameters()[Trk::theta], Trk::theta ) ;
+ std::vector<Trk::DefinedParameter> defPar ;
+ defPar.push_back( z0 ) ;
+ defPar.push_back( theta ) ;
+ if( !mp->covariance() ) return 0;
+ Trk::LocalParameters parFromSi( defPar ) ;
+ AmgSymMatrix(1) covFromSi;
+
+ covFromSi( 0, 0 ) = (*mp->covariance())( Trk::z0,Trk::z0 ) ;
+ covFromSi( 1, 1 ) = (*mp->covariance())( Trk::theta,Trk::theta ) ;
+ covFromSi( 1, 0 ) = (*mp->covariance())( Trk::z0, Trk::theta ) ;
+ covFromSi( 0, 1 ) = (*mp->covariance())( Trk::z0, Trk::theta ) ;
+
+ const Trk::Surface& mpSurf = mp->associatedSurface() ;
+
+ Trk::PseudoMeasurementOnTrack *pm = new Trk::PseudoMeasurementOnTrack( parFromSi
+ , covFromSi
+ , mpSurf
+ ) ;
+ return pm ;
+}
+
+const Trk::MeasurementSet RefitTracksAndVertex::addPM( Trk::MeasurementSet& ms, const Trk::PseudoMeasurementOnTrack* pm ) {
+ Trk::MeasurementSet sortedMS ;
+ sortedMS.push_back( pm ) ;
+ for( int i=0, i_max=ms.size() ; i!=i_max ; ++i ) {
+ sortedMS.push_back( ms[i] ) ;
+ }
+ return sortedMS ;
+}
+
+Trk::Track* RefitTracksAndVertex::fitSCTOnlyTrack(const Trk::Track* track) {
+ Trk::MeasurementSet setSCT ;
+ Trk::MeasurementSet setPix ;
+ Trk::MeasurementSet setTRT ;
+
+ const Trk::Perigee* perTrk = track->perigeeParameters();
+ if( !perTrk ) {
+ ATH_MSG_WARNING("RefitTracksAndVertex() : No Perigee parameter on track!");
+ return NULL ;
+ }
+
+ if( perTrk->eta() < m_selEtaMin || perTrk->eta() > m_selEtaMax)
+ return NULL ;
+
+ if( perTrk->pT() < m_selPtMin)
+ return NULL ;
+
+
+ //store all silicon measurements into the measurementset
+ DataVector<const Trk::MeasurementBase>::const_iterator it = track->measurementsOnTrack()->begin();
+ DataVector<const Trk::MeasurementBase>::const_iterator itEnd = track->measurementsOnTrack()->end();
+ for ( ; it!=itEnd; ++it){
+ if( !(*it) ) {
+ // log (MSG::WARNING) << "The MeasurementBase set has a void"
+ // << " member! Skip it.." << endreq;
+ } else {
+ const Trk::RIO_OnTrack* rio = dynamic_cast <const Trk::RIO_OnTrack*>(*it);
+ if (rio != 0) {
+ const Identifier& surfaceID = (rio->identify()) ;
+ if( m_idHelper->is_sct(surfaceID) ) {
+ setSCT.push_back ( *it ) ;
+ } else if( m_idHelper->is_trt(surfaceID) ) {
+ setTRT.push_back ( *it ) ;
+ } else if( m_idHelper->is_pixel(surfaceID) ){
+ setPix.push_back ( *it ) ;
+ }
+ }
+ }
+ }
+
+ if( (int)setSCT.size() < m_selNHitSCTMin )
+ return NULL;
+
+ ATH_MSG_DEBUG("RefitTracksAndVertex() : Found " << setSCT.size() << " SCT measurm's!" ) ;
+ ATH_MSG_DEBUG("RefitTracksAndVertex() : Found " << setPix.size() << " Pix measurm's!" ) ;
+ ATH_MSG_DEBUG("RefitTracksAndVertex() : Found " << setTRT.size() << " TRT measurm's!") ;
+
+
+ ATH_MSG_DEBUG( std::setiosflags( std::ios::scientific )) ;
+ ATH_MSG_DEBUG (std::setprecision( 7 )) ;
+ ATH_MSG_VERBOSE( std::setiosflags( std::ios::scientific )) ;
+ ATH_MSG_VERBOSE( std::setprecision( 7 )) ;
+
+ // now add z_0 and theta from original track as PseudoMeas to the TRT MeasurementSet
+ if(m_addPM){
+ const Trk::PseudoMeasurementOnTrack *pmFromSi = createPMfromSi( perTrk ) ;
+ if( !pmFromSi ) {
+ ATH_MSG_ERROR("RefitTracksAndVertex() : PseudoMeasurementOnTrack creation failed! " );
+ return NULL ;
+ }
+ ATH_MSG_DEBUG( "RefitTracksAndVertex() : pmFromSi " << *pmFromSi) ;
+ Trk::MeasurementSet setSCT = addPM( setSCT, pmFromSi ) ;
+ }
+
+ // Add TRT hits as they do no harm
+ for( int i=0, i_max=setTRT.size() ; i!=i_max ; ++i ) {
+ setSCT.push_back( setTRT[i] ) ;
+ }
+
+ ATH_MSG_VERBOSE ( "RefitTracksAndVertex() : Si+PM MeasurementSet : " );
+ for( int i=0, i_max=setSCT.size() ; i!=i_max ; ++i ) {
+ ATH_MSG_VERBOSE ("============== i=" << i << " =============");
+ ATH_MSG_VERBOSE ( *(setSCT[i]));
+ }
+ ATH_MSG_VERBOSE ("==========================================");
+
+ // fit TRT part of the track with PseudoMeas on z_0, theta
+ Trk::Track* trkSCT=m_trackFitter->fit( setSCT
+ , *perTrk
+ , true
+ , Trk::pion
+ //, Trk::muon
+ //, Trk::nonInteracting
+ ) ;
+ if( !trkSCT ) {
+ ATH_MSG_DEBUG( "RefitTracksAndVertex() : Fit of SCT part of the track failed! " ) ;
+ return NULL ;
+ }
+ const Trk::Perigee* perSCT = trkSCT->perigeeParameters();
+ ATH_MSG_VERBOSE( "RefitTracksAndVertex() : perSCT " << *perSCT) ;
+
+ return trkSCT;
+}
+
+MsgStream& operator<<( MsgStream& outst, const RefitTracksAndVertex& alg ) {
+ return alg.dump( outst ) ;
+}
+
+MsgStream& RefitTracksAndVertex::dump( MsgStream& outst ) const {
+ outst << std::endl ;
+ outst << "|-------------------------------------------------------------------";
+ outst << "-----------------------------|" << std::endl ;
+ outst << "| processed : "
+ << std::setw(7) << m_nTracksProcessed
+ << " tracks |"
+ << std::endl ;
+ outst << "| accepted by track presel. : "
+ << std::setw(7) << m_nTracksPresel
+ << " tracks |"
+ << std::endl ;
+ outst << "| accepted by track presel. + PM : "
+ << std::setw(7) << m_nTracksAccepted
+ << " tracks |"
+ << std::endl ;
+ outst << "| ------------------------------------------------------------------";
+ outst << "---------------------------- |" << std::endl ;
+ outst << "| reject by # PIX hits : "
+ << std::setw(7) << m_nRejectPIX
+ << " tracks |"
+ << std::endl ;
+ outst << "| reject by # SCT hits : "
+ << std::setw(7) << m_nRejectSCT
+ << " tracks |"
+ << std::endl ;
+ outst << "| reject by # TRT hits : "
+ << std::setw(7) << m_nRejectTRT
+ << " tracks |"
+ << std::endl ;
+ outst << "| ------------------------------------------------------------------";
+ outst << "---------------------------- |" << std::endl ;
+ outst << "| reject by exist. PM(TRT) : "
+ << std::setw(7) << m_nRejectPM
+ << " tracks |"
+ << std::endl ;
+ outst << "|-------------------------------------------------------------------";
+ outst << "-----------------------------|" << std::endl ;
+ return outst ;
+}
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.h b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.h
new file mode 100644
index 000000000000..fe45829816f2
--- /dev/null
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.h
@@ -0,0 +1,90 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// @file RefitTracksAndVertex.cxx
+//
+
+//
+// $Id: RefitTracksAndVertex.h,v 1.8 2009-05-18 15:48:05 moles Exp $
+
+#ifndef RefitTracksAndVertex_H
+#define RefitTracksAndVertex_H
+
+#include <string>
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "AthenaBaseComps/AthMsgStreamMacros.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "TrkFitterUtils/FitterTypes.h"
+#include "TrkParameters/TrackParameters.h"
+
+class AtlasDetectorID ;
+
+namespace Trk {
+ class ITrackFitter ;
+ class IVertexFitter ;
+ class Track ;
+ class TrackStateOnSurface ;
+ class PseudoMeasurementOnTrack ;
+}
+
+class RefitTracksAndVertex: public AthAlgorithm {
+ public:
+ /** The RefitTracksAndVertex is an implementation to add the so-called TRT momentum constraint on a track.
+ It was developed in the context of alignment, and its main idea is that if the TRT is free from "weak modes",
+ its curvature measurement can be added as a PseudoMeasurement (PM) to a Si-only track. This extra information,
+ taken into account while refitting the Si-only part of the track or read in by the alignment algorithms
+ directly, can serve to eliminate the alignment weak modes from the Si tracker of ATLAS.
+ More details can be found at: https://twiki.cern.ch/twiki/bin/view/Atlas/TRTMomConstraint */
+ RefitTracksAndVertex(const std::string& name, ISvcLocator* pSvcLocator) ;
+ ~RefitTracksAndVertex() ;
+ StatusCode initialize() ; //!< initialize method of this algorithm.
+ StatusCode execute() ; //!< execute method of this algorithm that is called for each event
+ StatusCode finalize() ; //!< finalize method of this algorithm. Prints out a summary of all events
+
+ MsgStream& dump( MsgStream& outst ) const ;
+
+ private:
+ // helper functions
+ /** Verifies if the given track passes the track selection criteria specified via the jobOptions */
+ bool accept( const Trk::Track& track ) ;
+ /** adds a PseudoMeasurement to a MeasurementSet */
+ const Trk::MeasurementSet addPM( Trk::MeasurementSet& ms, const Trk::PseudoMeasurementOnTrack* pm ) ;
+ /** creates a PseudoMeasurement with (z0, theta) from extended track perigee parameters */
+ const Trk::PseudoMeasurementOnTrack* createPMfromSi ( const Trk::Perigee* mp ) ;
+ /** Strips of all TRT hits from the track and replaces them with a TRT momentum constraint
+ from a the TRT segment belonging to the (extended) track (returns NULL in case of failure) */
+ Trk::Track* fitSCTOnlyTrack( const Trk::Track* track ) ;
+
+ // infrastructure
+ const AtlasDetectorID* m_idHelper ; //!< Detector ID helper
+ std::string m_vertexListInput ; //!< Name of the TrackCollection (input)
+
+
+
+ std::string m_trackListOutput ; //!< Name of the TrackCollection (Output)
+ std::string m_outputVertexContainerName ; //!< Name of vertex container
+ ToolHandle<Trk::ITrackFitter> m_trackFitter ; //!< The TrackFitter to refit the tracks (segment, momentum constraint)
+ ToolHandle<Trk::IVertexFitter> m_vertexFitter ; //!< The TrackFitter to refit the tracks (segment, momentum constraint)
+
+ // algorithm-internal members
+ bool m_applyTrkSel ; //!< apply a selection on tracks or not
+ bool m_refitTracks ; //!< refitTracks
+ bool m_addPM ; //!< apply a pseudo measurement based on the original track (theta,z0)
+ double m_selPtMin ; //!< minimal pT cut value for the TrackSelection
+ double m_selEtaMin ; //!< minimal eta cut value for the TrackSelection
+ double m_selEtaMax ; //!< maximal eta cut value for the TrackSelection
+ int m_selNHitPIXMin ; //!< minimal number of PIX hits cut value for the TrackSelection
+ int m_selNHitSCTMin ; //!< minimal number of SCT hits cut value for the TrackSelection
+ size_t m_nTracksProcessed ; //!< Counter for number of tracks processed
+ size_t m_nTracksPresel ; //!< Counter for number of tracks passing the preselection
+ size_t m_nTracksAccepted ; //!< Counter for number of tracks passing the preselection and with PM
+ size_t m_nRejectPIX ; //!< Counter for number of tracks failing the min number of PIX hits req.
+ size_t m_nRejectSCT ; //!< Counter for number of tracks failing the min number of SCT hits req.
+ size_t m_nRejectTRT ; //!< Counter for number of tracks failing the min number of TRT hits req.
+ size_t m_nRejectPM ; //!< Counter for number of tracks failing the addition of a pseudo-measurement (PM)
+} ;
+
+MsgStream& operator<<( MsgStream& outst, const RefitTracksAndVertex& alg ) ;
+
+#endif
diff --git a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/components/InDetBeamSpotFinder_entries.cxx b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/components/InDetBeamSpotFinder_entries.cxx
index 7a155da2f5ce..66d496905bac 100644
--- a/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/components/InDetBeamSpotFinder_entries.cxx
+++ b/InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/components/InDetBeamSpotFinder_entries.cxx
@@ -1,14 +1,19 @@
#include "../InDetBeamSpotFinder.h"
#include "../InDetBeamSpotRooFit.h"
+#include "../InDetBeamSpotVertex.h"
+#include "../RefitTracksAndVertex.h"
+
#include "GaudiKernel/DeclareFactoryEntries.h"
DECLARE_NAMESPACE_ALGORITHM_FACTORY(InDet, InDetBeamSpotFinder )
DECLARE_NAMESPACE_TOOL_FACTORY(InDet, InDetBeamSpotVertex)
DECLARE_NAMESPACE_TOOL_FACTORY(InDet, InDetBeamSpotRooFit)
+DECLARE_ALGORITHM_FACTORY( RefitTracksAndVertex )
DECLARE_FACTORY_ENTRIES(InDetBeamSpotFinder){
DECLARE_NAMESPACE_ALGORITHM(InDet, InDetBeamSpotFinder)
DECLARE_NAMESPACE_TOOL(InDet, InDetBeamSpotVertex)
DECLARE_NAMESPACE_TOOL(InDet, InDetBeamSpotRooFit)
+ DECLARE_ALGORITHM(RefitTracksAndVertex)
}
--
GitLab