From 1d24bf713d2f7ea9c92b0fae8c1cdeeff0da6b45 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 31 Oct 2023 11:54:25 +0000
Subject: [PATCH] DecodedDataFromMCRichHits: Add 0-25 time window, option to
reject scintillation
---
.../component/DecodedDataFromMCRichHits.cpp | 92 ++++++++++++-------
1 file changed, 57 insertions(+), 35 deletions(-)
diff --git a/Rich/RichFutureMCUtils/src/component/DecodedDataFromMCRichHits.cpp b/Rich/RichFutureMCUtils/src/component/DecodedDataFromMCRichHits.cpp
index d8067b837c8..ac7ce7f6c4c 100644
--- a/Rich/RichFutureMCUtils/src/component/DecodedDataFromMCRichHits.cpp
+++ b/Rich/RichFutureMCUtils/src/component/DecodedDataFromMCRichHits.cpp
@@ -9,6 +9,10 @@
* or submit itself to any jurisdiction. *
\*****************************************************************************/
+// Gaudi (must be first)
+#include "GaudiKernel/ParsersFactory.h"
+#include "GaudiKernel/StdArrayAsProperty.h"
+
// Rich (Future) Kernel
#include "RichFutureKernel/RichAlgBase.h"
@@ -68,7 +72,8 @@ namespace Rich::Future::MC {
auto& summaries = std::get<LHCb::MCRichDigitSummarys>( out_data );
// map from Smart ID to associated MC hits
- std::unordered_map<LHCb::RichSmartID, std::vector<LHCb::MCRichHit*>> hitsPerPix;
+ using HitData = std::pair<const LHCb::MCRichHit*, double>;
+ std::unordered_map<LHCb::RichSmartID, std::vector<HitData>> hitsPerPix;
// Get the location of the container of an MC hit in the TES
auto hitLocation = [&]( const auto* obj ) { return ( obj ? objectLocation( obj->parent() ) : "Not Contained" ); };
@@ -77,20 +82,16 @@ namespace Rich::Future::MC {
auto getTime = [&]( const auto* mchit ) {
if ( !mchit ) { return 0.0; }
// Deduce if this is a spillover event and if it is apply offset
- const auto loc = hitLocation( mchit );
- double spillOffset = 0.0;
- if ( loc.find( "PrevPrev" ) != std::string::npos ) {
- spillOffset = -50.0;
- } else if ( loc.find( "Prev" ) != std::string::npos ) {
- spillOffset = -25.0;
- } else if ( loc.find( "NextNext" ) != std::string::npos ) {
- spillOffset = 50.0;
- } else if ( loc.find( "Next" ) != std::string::npos ) {
- spillOffset = 25.0;
- }
- _ri_verbo << "Loc " << loc << " " << spillOffset << endmsg;
- // return offset + TOF
- return spillOffset + mchit->timeOfFlight();
+ const auto loc = hitLocation( mchit );
+ const double spillOffset = //
+ ( ( loc.find( "PrevPrev" ) != std::string::npos ) ? -50.0 : //
+ ( loc.find( "Prev" ) != std::string::npos ) ? -25.0 : //
+ ( loc.find( "NextNext" ) != std::string::npos ) ? 50.0 : //
+ ( loc.find( "Next" ) != std::string::npos ) ? 25.0 : //
+ 0.0 );
+ const auto hitT = spillOffset + mchit->timeOfFlight();
+ _ri_verbo << "TES Location '" << loc << "' SpillOffset=" << spillOffset << " HitTime=" << hitT << endmsg;
+ return hitT;
};
// process a container of MC hits
@@ -98,14 +99,27 @@ namespace Rich::Future::MC {
for ( const auto mchit : mchits ) {
const auto id = mchit->sensDetID();
if ( id.isValid() ) {
+ _ri_verbo << id << " " << mchit->mcRichDigitHistoryCode() << endmsg;
// Apply all sorts MC based filtering here, e.g. signal only, timing windows etc.
// To be extended as needed
- if ( m_rejectBackground && mchit->isBackground() ) {
- _ri_verbo << id << " Rejected as background" << endmsg;
+ if ( m_rejectAllBackgrounds && mchit->isBackground() ) {
+ _ri_verbo << " -> Rejected as background" << endmsg;
continue;
}
- // finally save this id -> hit mapping
- hitsPerPix[id].push_back( mchit );
+ if ( m_rejectScintillation && mchit->radScintillation() ) {
+ _ri_verbo << " -> Rejected as scintillation" << endmsg;
+ continue;
+ }
+ // Are we inside the 0-25 ns window ?
+ const auto hitT = getTime( mchit );
+ // time with shift for each RICH
+ const auto hitTShifted = hitT - m_timeShift[id.rich()];
+ _ri_verbo << "Shifted time = " << hitTShifted << endmsg;
+ if ( hitTShifted >= 0.0 && hitTShifted <= 25.0 ) {
+ _ri_verbo << " -> In Time" << endmsg;
+ // finally save this id -> hit mapping
+ hitsPerPix[id].emplace_back( mchit, hitT );
+ }
} // hit is valid
} // mc hit loop
};
@@ -131,6 +145,9 @@ namespace Rich::Future::MC {
summaries.reserve( hitsPerPix.size() );
for ( const auto& [id, hits] : hitsPerPix ) {
+ // Should never fail but just in case...
+ if ( hits.empty() ) { continue; }
+
// temporary copy of SmartID
auto newID = id;
@@ -151,38 +168,36 @@ namespace Rich::Future::MC {
// Add best time for this ID
if ( m_includeTimeInfo ) {
// attempt to select the 'best' time for this ID ...
- const LHCb::MCRichHit* mch = nullptr;
- if ( !hits.empty() ) { // should never fail ...
- mch = hits.front();
- for ( const auto hit : hits ) {
- if ( hit->isSignal() ) {
- mch = hit;
- break;
- }
+ const auto* mch = &hits.front();
+ for ( const auto& hit : hits ) {
+ if ( hit.first->isSignal() ) {
+ mch = &hit;
+ break;
}
}
- newID.setTime( getTime( mch ) );
+ newID.setTime( mch->second );
}
// Add (dd4hep corrected) ID to decoded list
ids.emplace_back( newID );
// Form MC summaries for each MCHit
- for ( const auto hit : hits ) {
+ for ( const auto& hit : hits ) {
auto summary = std::make_unique<LHCb::MCRichDigitSummary>();
// do not want dd4hep corrected ID here as utility does its own correction.
summary->setRichSmartID( id );
- summary->setMCParticle( hit->mcParticle() );
- auto history = hit->mcRichDigitHistoryCode();
+ summary->setMCParticle( hit.first->mcParticle() );
+ auto history = hit.first->mcRichDigitHistoryCode();
history.setSignalEvent( true ); // normally done by Boole for signal event only
summary->setHistory( history );
summaries.add( summary.release() );
}
- }
+
+ } // loop over hits
// Sort
SmartIDSorter::sortByRegion( ids );
- for ( const auto id : ids ) { ri_debug( id, endmsg ); }
+ // for ( const auto id : ids ) { ri_debug( id, endmsg ); }
// fill decoded data structure
decoded_data.addSmartIDs( ids );
@@ -213,11 +228,18 @@ namespace Rich::Future::MC {
/// Include time information in the decoded data
Gaudi::Property<bool> m_includeTimeInfo{this, "IncludeTimeInfo", true};
- /// Reject background hits
- Gaudi::Property<bool> m_rejectBackground{this, "RejectBackground", false};
+ /// Reject all background hits
+ Gaudi::Property<bool> m_rejectAllBackgrounds{this, "RejectAllBackgrounds", false};
+
+ /// Reject scintillation hits
+ Gaudi::Property<bool> m_rejectScintillation{this, "RejectScintillation", false};
/// Include spillover events
Gaudi::Property<bool> m_enableSpillover{this, "EnableSpillover", true};
+
+ /// Time window shift for each RICH
+ Gaudi::Property<DetectorArray<double>> m_timeShift{
+ this, "TimeCalib", {0., 40.}, "Global time shift for each RICH, to get both to same calibrated point"};
};
// Declaration of the Algorithm Factory
--
GitLab