diff --git a/Event/xAOD/xAODEventInfo/Root/EventInfo_v1.cxx b/Event/xAOD/xAODEventInfo/Root/EventInfo_v1.cxx
index 4ce5ff397c76d39b7d0397616fba8151c7768f6a..43a8c8e2ab7e771aa994145910500da15a14d211 100644
--- a/Event/xAOD/xAODEventInfo/Root/EventInfo_v1.cxx
+++ b/Event/xAOD/xAODEventInfo/Root/EventInfo_v1.cxx
@@ -949,6 +949,27 @@ namespace xAOD {
AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( EventInfo_v1, uint32_t, beamStatus,
setBeamStatus )
+
+ /// Accessor for "BeamSpotWeight"
+ static const SG::AuxElement::Accessor< float >
+ BeamSpotWeight( "BeamSpotWeight" );
+
+ bool EventInfo_v1::hasBeamSpotWeight() const {
+ return BeamSpotWeight.isAvailable( *this );
+ }
+
+ float EventInfo_v1::beamSpotWeight() const {
+ // If the value is not available, then return 1.0
+ if( ! BeamSpotWeight.isAvailable( *this ) ) return 1.0f;
+ return BeamSpotWeight( *this );
+ }
+
+ void EventInfo_v1::setBeamSpotWeight( float value ) {
+ BeamSpotWeight( *this ) = value;
+ return;
+ }
+
+
//
/////////////////////////////////////////////////////////////////////////////
diff --git a/Event/xAOD/xAODEventInfo/xAODEventInfo/versions/EventInfo_v1.h b/Event/xAOD/xAODEventInfo/xAODEventInfo/versions/EventInfo_v1.h
index dbc7ee3fa74304ae1eaa36ee130ec11b2d30f183..88e334702ebe58bcd26944b83a44f36aec637f29 100644
--- a/Event/xAOD/xAODEventInfo/xAODEventInfo/versions/EventInfo_v1.h
+++ b/Event/xAOD/xAODEventInfo/xAODEventInfo/versions/EventInfo_v1.h
@@ -1,7 +1,7 @@
// Dear emacs, this is -*- c++ -*-
/*
- Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
// $Id: EventInfo_v1.h 727083 2016-03-01 15:20:50Z krasznaa $
@@ -484,6 +484,13 @@ namespace xAOD {
/// Set the beam spot's status word
void setBeamStatus( uint32_t value );
+ /// Weight for beam spot size reweighting
+ float beamSpotWeight() const;
+ /// Check if weight for beam spot size reweighting exists
+ bool hasBeamSpotWeight() const;
+ /// Set weight for beam spot size reweighting
+ void setBeamSpotWeight( float value );
+
/// @}
/// @name Functions used by pile-up digitisation
diff --git a/Simulation/BeamEffects/src/BeamSpotReweightingAlg.cxx b/Simulation/BeamEffects/src/BeamSpotReweightingAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fc8a8c89f806f496d8fa7985a3e6400ccbbc665d
--- /dev/null
+++ b/Simulation/BeamEffects/src/BeamSpotReweightingAlg.cxx
@@ -0,0 +1,180 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+// class header
+#include "BeamSpotReweightingAlg.h"
+
+// Athena headers
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/WriteHandle.h"
+#include "StoreGate/WriteDecorHandle.h"
+
+// HepMC includes
+#include "AtlasHepMC/GenEvent.h"
+#include "AtlasHepMC/GenVertex.h"
+
+// CLHEP includes
+#include "CLHEP/Vector/LorentzVector.h"
+
+namespace Simulation
+{
+
+ BeamSpotReweightingAlg::BeamSpotReweightingAlg( const std::string& name, ISvcLocator* pSvcLocator )
+ : AthReentrantAlgorithm( name, pSvcLocator )
+ {
+ declareProperty( "Input_beam_sigma_z", m_input_beam_sigma_z, "Beam spot sigma of the input HIT file to be reweighted");
+ }
+
+ /** Athena algorithm's interface method initialize() */
+ StatusCode BeamSpotReweightingAlg::initialize()
+ {
+ // This will check that the properties were initialized
+ // properly by job configuration.
+ ATH_CHECK( m_inputMcEventCollection.initialize() );
+ ATH_CHECK( m_beamSpotWeight.initialize() );
+ ATH_CHECK( m_beamSpotKey.initialize() );
+ ATH_CHECK( m_eventInfo.initialize() );
+
+ ATH_MSG_INFO("Input_beam_sigma_z="<<m_input_beam_sigma_z);
+
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode BeamSpotReweightingAlg::execute(const EventContext& ctx) const
+ {
+ float weight=1;
+ if(m_input_beam_sigma_z<0) {
+ ATH_MSG_DEBUG("Input_beam_sigma_z="<<m_input_beam_sigma_z<<"<0, do nothing");
+ return StatusCode::SUCCESS;
+ }
+
+ if(msgLvl(MSG::DEBUG)) {
+ SG::ReadHandle<xAOD::EventInfo> evt (m_eventInfo,ctx);
+ if (!evt.isValid()) {
+ ATH_MSG_FATAL( "Could not find event info" );
+ return(StatusCode::FAILURE);
+ }
+ ATH_MSG_DEBUG( "EventInfo before adding the weight: " << evt->eventNumber()
+ << " run: " << evt->runNumber()
+ << " hasBeamSpotWeight: "<< evt->hasBeamSpotWeight()
+ << " beamSpotWeight: "<< evt->beamSpotWeight() );
+ }
+
+ // Construct handles from the keys.
+ SG::ReadHandle<McEventCollection> h_inputMcEventCollection (m_inputMcEventCollection, ctx);
+ if(!h_inputMcEventCollection.isValid()) {
+ ATH_MSG_FATAL("No input McEventCollection called " << h_inputMcEventCollection.name() << " in StoreGate.");
+ return StatusCode::FAILURE;
+ }
+
+ // loop over the event in the mc collection
+ for (const auto& currentGenEvent : *h_inputMcEventCollection) {
+ // skip empty events
+ if ( !currentGenEvent ) {
+ //the hard scatter event is always the first one. If not present, do nothing
+ ATH_MSG_WARNING("No first event found in the McEventCollection input collection, use default weight="<<weight);
+ break;
+ }
+ HepMC::GenVertex *signalVertex = GetSignalProcessVertex(*currentGenEvent);
+
+ if(signalVertex) {
+ //now calculate weight
+ weight=1.0;
+ float z=signalVertex->position().z();
+ SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
+ float beamz=beamSpotHandle->beamPos().z();
+ float newsigmaz=beamSpotHandle->beamSigma(2);
+ float pullold=(z-beamz)/m_input_beam_sigma_z;
+ float pullnew=(z-beamz)/newsigmaz;
+
+ ATH_MSG_DEBUG("Beamspot z="<<beamz<<", z="<<z<<"; old sigma_z="<<m_input_beam_sigma_z<<", old pull="<<pullold<<"; new sigma_z="<<newsigmaz<<", new pull="<<pullnew);
+
+ if(std::abs(pullold)<10 && std::abs(pullnew)<10) {
+ //Use Gauss probability ratio to calculate the weight:
+ //weight=TMath::Gaus(z,0,35,true)/TMath::Gaus(z,0,42,true);
+ //This can be simplified to:
+ weight=m_input_beam_sigma_z/newsigmaz * std::exp(0.5*(pullold*pullold - pullnew*pullnew));
+ } else {
+ ATH_MSG_WARNING("Large pull of beamspot: Beamspot z="<<beamz<<", z="<<z<<"; old sigma_z="<<m_input_beam_sigma_z<<", old pull="<<pullold<<"; new sigma_z="<<newsigmaz<<", new pull="<<pullnew<<" => use default weight="<<weight);
+ }
+ } else {
+ ATH_MSG_WARNING("No signal vertex found, use default weight="<<weight);
+ }
+
+ ATH_MSG_DEBUG("Beam weight="<<weight);
+ break;
+ }
+
+ SG::WriteDecorHandle<xAOD::EventInfo,float> BeamSpotWeight(m_beamSpotWeight,ctx);
+ if (!BeamSpotWeight.isPresent()) {
+ ATH_MSG_ERROR( "BeamSpotWeight.isPresent check fails" );
+ return StatusCode::FAILURE;
+ }
+ if (!BeamSpotWeight.isAvailable()) {
+ BeamSpotWeight(0) = weight;
+ }
+
+ if(msgLvl(MSG::DEBUG)) {
+ SG::ReadHandle<xAOD::EventInfo> evt (m_eventInfo,ctx);
+ if (!evt.isValid()) {
+ ATH_MSG_FATAL( "Could not find event info" );
+ return(StatusCode::FAILURE);
+ }
+ ATH_MSG_DEBUG( "EventInfo after adding the weight: " << evt->eventNumber()
+ << " run: " << evt->runNumber()
+ << " hasBeamSpotWeight: "<< evt->hasBeamSpotWeight()
+ << " beamSpotWeight: "<< evt->beamSpotWeight() );
+ }
+
+ return StatusCode::SUCCESS;
+ }
+
+ HepMC::GenVertex* BeamSpotReweightingAlg::GetSignalProcessVertex(const HepMC::GenEvent& ge) const
+ {
+ //Ensure that we have a valid signal_process_vertex
+#ifdef HEPMC3
+ if( !HepMC::signal_process_vertex(ge) ) {
+ if (!ge.vertices_empty()) {
+ ATH_MSG_DEBUG("No signal_process_vertex found - using the first GenVertex in the event.");
+ HepMC::GenVertexPtr signalVertex = ge.vertices().front();
+ return signalVertex;
+ }
+ if( !HepMC::signal_process_vertex(ge) ) { // Insanity check
+ if (!ge.vertices().empty()) {
+ ATH_MSG_ERROR("Failed to set signal_process_vertex for GenEvent!!");
+ return nullptr;
+ }
+ ATH_MSG_WARNING("No signal_process_vertex found. Empty GenEvent!");
+ return nullptr;
+ }
+ }
+ else {
+ ATH_MSG_DEBUG("signal_process_vertex set by Generator.");
+ return HepMC::signal_process_vertex(ge);
+ }
+#else
+ if( !ge.signal_process_vertex() ) {
+ if (!ge.vertices_empty()) {
+ ATH_MSG_DEBUG("No signal_process_vertex found - using the first GenVertex in the event.");
+ HepMC::GenVertex *signalVertex = *(ge.vertices_begin());
+ return signalVertex;
+ }
+ if( !ge.signal_process_vertex() ) { // Insanity check
+ if (!ge.vertices_empty()) {
+ ATH_MSG_ERROR("Failed to set signal_process_vertex for GenEvent!!");
+ return nullptr;
+ }
+ ATH_MSG_WARNING("No signal_process_vertex found. Empty GenEvent!");
+ return nullptr;
+ }
+ }
+ else {
+ ATH_MSG_DEBUG("signal_process_vertex set by Generator.");
+ return ge.signal_process_vertex();
+ }
+#endif
+ return nullptr;
+ }
+
+} // namespace Simulation
diff --git a/Simulation/BeamEffects/src/BeamSpotReweightingAlg.h b/Simulation/BeamEffects/src/BeamSpotReweightingAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..83c1f4141be18cf7cf2c8c186261f6477cddd6db
--- /dev/null
+++ b/Simulation/BeamEffects/src/BeamSpotReweightingAlg.h
@@ -0,0 +1,75 @@
+// Dear emacs, this is -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef BEAMEFFECTS_BeamSpotReweightingAlg_H
+#define BEAMEFFECTS_BeamSpotReweightingAlg_H
+
+// base class header
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Athena includes
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+#include "StoreGate/WriteDecorHandleKey.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "BeamSpotConditionsData/BeamSpotData.h"
+#include "GeneratorObjects/McEventCollection.h"
+
+// Gaudi includes
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/PhysicalConstants.h"
+
+// Forward declarations
+#include "AtlasHepMC/GenEvent_fwd.h"
+
+namespace HepMC
+{
+ class GenVertex;
+}
+
+namespace Simulation
+{
+
+ /** @class BeamSpotReweightingAlg
+
+ This AthReentrantAlgorithm calculates an event weight based such that
+ the distribution of z-positions of events matches a new beamspot size
+ */
+ class BeamSpotReweightingAlg : public AthReentrantAlgorithm
+ {
+ public:
+
+ //** Constructor with parameters */
+ BeamSpotReweightingAlg( const std::string& name, ISvcLocator* pSvcLocator );
+
+ /** Destructor */
+ virtual ~BeamSpotReweightingAlg() = default;
+
+ /** Athena algorithm's interface method initialize() */
+ virtual StatusCode initialize() override final;
+
+ /** Athena algorithm's interface method execute() */
+ virtual StatusCode execute(const EventContext& ctx) const override final;
+
+ private:
+
+ /** Ensure that the GenEvent::signal_process_vertex has been set */
+ HepMC::GenVertex* GetSignalProcessVertex(const HepMC::GenEvent& ge) const;
+
+ SG::ReadHandleKey<McEventCollection> m_inputMcEventCollection { this, "InputMcEventCollection", "TruthEvent", "The name of the input McEventCollection" };
+
+ SG::WriteDecorHandleKey<xAOD::EventInfo> m_beamSpotWeight {this, "BeamSpotWeight", "EventInfo.BeamSpotWeight", "Decoration for a beam spot weight when reweighting the beam spot size" };
+
+ SG::ReadCondHandleKey<InDet::BeamSpotData> m_beamSpotKey { this, "BeamSpotKey", "BeamSpotData", "SG key for beam spot" };
+
+ SG::ReadHandleKey<xAOD::EventInfo> m_eventInfo { this, "EventInfo", "EventInfo", "EventInfo object in SG" };
+
+ float m_input_beam_sigma_z = 42*Gaudi::Units::mm;
+ };
+
+}
+
+#endif // BEAMEFFECTS_BeamSpotReweightingAlg_H
diff --git a/Simulation/BeamEffects/src/components/BeamEffects_entries.cxx b/Simulation/BeamEffects/src/components/BeamEffects_entries.cxx
index 019e319d30a7fe193567daefc1d21d6b14e951ea..9ed156c684dc2b0e5d785aee74e8470f1aeac360 100644
--- a/Simulation/BeamEffects/src/components/BeamEffects_entries.cxx
+++ b/Simulation/BeamEffects/src/components/BeamEffects_entries.cxx
@@ -8,6 +8,7 @@
#include "../GenEventBeamEffectBooster.h"
#include "../GenEventRotator.h"
#include "../BeamEffectsAlg.h"
+#include "../BeamSpotReweightingAlg.h"
DECLARE_COMPONENT( Simulation::ZeroLifetimePositioner )
DECLARE_COMPONENT( Simulation::GenEventValidityChecker )
DECLARE_COMPONENT( Simulation::GenEventVertexPositioner )
@@ -18,4 +19,5 @@ DECLARE_COMPONENT( Simulation::VertexPositionFromFile )
DECLARE_COMPONENT( Simulation::GenEventBeamEffectBooster )
DECLARE_COMPONENT( Simulation::GenEventRotator )
DECLARE_COMPONENT( Simulation::BeamEffectsAlg )
+DECLARE_COMPONENT( Simulation::BeamSpotReweightingAlg )