Updates for release 21 (InDetBeamSpotFinder-01-02-00)

  *  Update to allow:
	  * Fits done per unit time
		* Fix additional parameters
		* plb fits

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 )

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;
   };
 }

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

InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/BeamspotRefitVertex.py

+
+# 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
+

InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/share/postInclude.addSCTonlyToAODFile.py

+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.' )
+
+

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;

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

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){    

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 
 

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;
 }

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){

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
     

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;
 

InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.cxx

+/*
+  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 ;
+}

InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/RefitTracksAndVertex.h

+/*
+  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

InnerDetector/InDetCalibAlgs/InDetBeamSpotFinder/src/components/InDetBeamSpotFinder_entries.cxx

 #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)
 }