Commit d37d422f authored by Chao Peng's avatar Chao Peng
Browse files

General PMT Hit Digitization and Reconstruction

parent c2eba219
/* General PhotoMultiplier Digitization
*
* Apply the given quantum efficiency for photon detection
* Converts the number of detected photons to signal amplitude
*
* Author: Chao Peng (ANL)
* Date: 10/02/2020
*/
#include <iterator>
#include <algorithm>
#include <unordered_map>
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "JugBase/DataHandle.h"
// Event Model related classes
#include "eicd/RawPMTHitCollection.h"
#include "dd4pod/PhotoMultiplierHitCollection.h"
using namespace Gaudi::Units;
namespace Jug::Digi {
class PhotoMultiplierDigi : public GaudiAlgorithm
{
public:
DataHandle<dd4pod::PhotoMultiplierHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawPMTHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
Gaudi::Property<std::vector<std::pair<double, double>>>
u_quantumEfficiency{this, "quantumEfficiency", {{2.6*eV, 0.3}, {7.0*eV, 0.3}}};
Gaudi::Property<double> m_hitTimeWindow{this, "hitTimeWindow", 20.0*ns};
Gaudi::Property<double> m_timeStep{this, "timeStep", 0.0625*ns};
Gaudi::Property<double> m_speMean{this, "speMean", 80.0};
Gaudi::Property<double> m_speError{this, "speError", 16.0};
Gaudi::Property<double> m_pedMean{this, "pedMean", 200.0};
Gaudi::Property<double> m_pedError{this, "pedError", 3.0};
Rndm::Numbers m_rngUni, m_rngNorm;
// constructor
PhotoMultiplierDigi(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection,"");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
auto sc1 = m_rngUni.initialize(randSvc, Rndm::Flat(0., 1.));
auto sc2 = m_rngNorm.initialize(randSvc, Rndm::Gauss(0., 1.));
if (!sc1.isSuccess() || !sc2.isSuccess()) {
error() << "Cannot initialize random generator!" << endmsg;
return StatusCode::FAILURE;
}
qe_init();
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collection
const auto &sim = *m_inputHitCollection.get();
// Create output collections
auto &raw = *m_outputHitCollection.createAndPut();
struct HitData { int npe; double signal; double time; };
std::unordered_map<long long, std::vector<HitData>> hit_groups;
// collect the photon hit in the same cell
// calculate signal
for(const auto& ahit : sim) {
// quantum efficiency
if (!qe_pass(ahit.energy(), m_rngUni())) {
continue;
}
// cell id, time, signal amplitude
long long id = ahit.cellID();
double time = ahit.truth().time;
double amp = m_speMean + m_rngNorm()*m_speError;
// group hits
auto it = hit_groups.find(id);
if (it != hit_groups.end()) {
size_t i = 0;
for (auto git = it->second.begin(); git != it->second.end(); ++git, ++i) {
if (std::abs(time - git->time) <= m_hitTimeWindow) {
git->npe += 1;
git->signal += amp;
break;
}
}
// no hits group found
if (i >= it->second.size()) {
it->second.emplace_back(HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time});
}
} else {
hit_groups[id] = {HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}};
}
}
// build hit
for (auto &it : hit_groups) {
for (auto &data : it.second) {
eic::RawPMTHit hit(it.first, (unsigned) data.signal, (unsigned) data.time/m_timeStep);
raw.push_back(hit);
}
}
return StatusCode::SUCCESS;
}
private:
void qe_init()
{
auto &qeff = u_quantumEfficiency.value();
// sort quantum efficiency data first
std::sort(qeff.begin(), qeff.end(),
[] (const std::pair<double, double> &v1, const std::pair<double, double> &v2) {
return v1.first < v2.first;
});
// sanity checks
if (qeff.empty()) {
qeff = {{2.6*eV, 0.3}, {7.0*eV, 0.3}};
warning() << "Invalid quantum efficiency data provided, using default values: " << qeff << endmsg;
}
if (qeff.front().first > 3.0*eV) {
warning() << "Quantum efficiency data start from " << qeff.front().first/eV
<< " eV, maybe you are using wrong units?" << endmsg;
}
if (qeff.back().first < 6.0*eV) {
warning() << "Quantum efficiency data end at " << qeff.back().first/eV
<< " eV, maybe you are using wrong units?" << endmsg;
}
}
// helper function for linear interpolation
// Comp return is defined as: equal, 0; greater, > 0; less, < 0
template<class RndmIter, typename T, class Compare>
RndmIter interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp)
{
// special cases
auto dist = std::distance(beg, end);
if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) {
return end;
}
auto first = beg, last = end, mid = std::next(beg, dist/2);
while (mid != end) {
if (comp(*mid, val) == 0) {
return mid;
} else if (comp(*mid, val) > 0) {
end = mid;
} else {
beg = std::next(mid);
}
mid = std::next(beg, std::distance(beg, end)/2);
}
if (mid == last || comp(*mid, val) > 0) {
return std::prev(mid);
}
return mid;
}
bool qe_pass(double ev, double rand)
{
auto &qeff = u_quantumEfficiency.value();
auto it = interval_search(qeff.begin(), qeff.end(), ev,
[] (const std::pair<double, double> &vals, double val) {
return vals.first - val;
});
if (it == qeff.end()) {
// info() << ev/eV << " eV is out of QE data range, assuming 0% efficiency" << endmsg;
return false;
}
double prob = it->second;
auto itn = std::next(it);
if (itn != qeff.end() && (itn->first - it->first != 0)) {
prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first);
}
// info() << ev/eV << " eV, QE: " << prob*100. << "%" << endmsg;
return rand <= prob;
}
};
DECLARE_COMPONENT(PhotoMultiplierDigi)
} // namespace Jug::Digi
from Gaudi.Configuration import *
from GaudiKernel import SystemOfUnits as units
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput
from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
podioevent = EICDataSvc("EventDataSvc", inputs=["derp.root"], OutputLevel=DEBUG)
geo_service = GeoSvc("GeoSvc")
podioevent = EICDataSvc("EventDataSvc", inputs=["rich_test.root"], OutputLevel=DEBUG)
from Configurables import PodioInput
from Configurables import Jug__Digi__ExampleCaloDigi as ExampleCaloDigi
from Configurables import Jug__Digi__UFSDTrackerDigi as UFSDTrackerDigi
podioinput = PodioInput("PodioReader", collections=["mcparticles","LAEC_PrShHits","LAEC_ShHits","FAEC_PrShHits","FAEC_ShHits","GEMTrackerHits"], OutputLevel=DEBUG)
caldigi = ExampleCaloDigi(inputHitCollection="FAEC_ShHits",outputHitCollection="RawFAECShowerHits")
ufsd_digi = UFSDTrackerDigi(inputHitCollection="GEMTrackerHits",outputHitCollection="GEMRawHits")
from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
out = PodioOutput("out", filename="test.root")
qe_data = [(1.0, 0.25), (7.5, 0.25),]
podioinput = PodioInput("PodioReader", collections=["mcparticles", "ForwardRICHHits"], OutputLevel=DEBUG)
pmtdigi = PhotoMultiplierDigi(inputHitCollection="ForwardRICHHits", outputHitCollection="DigiForwardRICHHits",
quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
out = PodioOutput("out", filename="digi_rich_test.root")
out.outputCommands = ["keep *"]
ApplicationMgr(
TopAlg = [podioinput, caldigi,ufsd_digi, out
],
TopAlg = [podioinput, pmtdigi, out],
EvtSel = 'NONE',
EvtMax = 5,
EvtMax = 100,
ExtSvc = [podioevent],
OutputLevel=DEBUG
)
......
......@@ -37,6 +37,7 @@ gaudi_add_module(JugRecoPlugins
src/components/CalorimeterIslandCluster.cpp
src/components/CrystalEndcapsReco.cpp
src/components/ClusterRecoCoG.cpp
src/components/PhotoMultiplierReco.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
......
/* General PhotoMultiplier Reconstruction
*
* Apply the given quantum efficiency for photon detection
* Converts the number of s to signal amplitude
*
* Author: Chao Peng (ANL)
* Date: 10/03/2020
*/
#include <algorithm>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
// FCCSW
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
// Event Model related classes
#include "eicd/PMTHitCollection.h"
#include "eicd/RawPMTHitCollection.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
class PhotoMultiplierReco : public GaudiAlgorithm
{
public:
DataHandle<eic::RawPMTHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::PMTHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
Gaudi::Property<double> m_timeStep{this, "timeStep", 0.0625*ns};
Gaudi::Property<double> m_minNpe{this, "minNpe", 0.5};
Gaudi::Property<double> m_speMean{this, "speMean", 80.0};
Gaudi::Property<double> m_pedMean{this, "pedMean", 200.0};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
PhotoMultiplierReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service("GeoSvc");
if (!m_geoSvc) {
error() << "Unable to locate Geometry Service. "
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto &rawhits = *m_inputHitCollection.get();
// Create output collections
auto &hits = *m_outputHitCollection.createAndPut();
// reconstrut number of photo-electrons and time
for (auto &rh : rawhits) {
float npe = (rh.amplitude() - m_pedMean)/m_speMean;
if (npe >= m_minNpe) {
float time = rh.timeStamp()*m_timeStep;
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
hits.push_back(eic::PMTHit{
id, npe, time,
{gpos.x(), gpos.y(), gpos.z()},
{pos.x(), pos.y(), pos.z()}
});
}
}
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(PhotoMultiplierReco)
} // namespace Jug::Reco
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment