diff --git a/generators/Herwig/HepMCTrigger/HepMCJetTrigger.cc b/generators/Herwig/HepMCTrigger/HepMCJetTrigger.cc index 02ca208f25..792720a962 100644 --- a/generators/Herwig/HepMCTrigger/HepMCJetTrigger.cc +++ b/generators/Herwig/HepMCTrigger/HepMCJetTrigger.cc @@ -96,6 +96,11 @@ std::vector HepMCJetTrigger::findAllJets(HepMC::GenEvent* e1 if (!(*iter)->end_vertex() && (*iter)->status() == 1) { auto p = (*iter)->momentum(); + auto pd = std::abs((*iter)->pdg_id()); + if (pd >= 12 && pd <= 18) + { + continue; // keep jet in the expected behavioro + } fastjet::PseudoJet pj(p.px(), p.py(), p.pz(), p.e()); pj.set_user_index((*iter)->barcode()); input.push_back(pj); @@ -122,6 +127,10 @@ int HepMCJetTrigger::jetsAboveThreshold(const std::vector& j for (const auto& j : jets) { float const pt = j.pt(); + if (std::abs(j.eta()) > 1.1) + { + continue; + } if (pt > this->threshold) { n_good_jets++; diff --git a/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.cc b/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.cc new file mode 100644 index 0000000000..fe1193cba0 --- /dev/null +++ b/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.cc @@ -0,0 +1,369 @@ +#include "HepMCParticleTrigger.h" + +#include +#include +#include + +#include + +#include +#include + +#include + +#include + +#include +#include +#include +#include +//____________________________________________________________________________.. +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +HepMCParticleTrigger::HepMCParticleTrigger(float trigger_thresh, int n_incom, bool up_lim, const std::string& name) + : SubsysReco(name) + , threshold(trigger_thresh) + , goal_event_number(n_incom) + , set_event_limit(up_lim) + , _theEtaHigh(1.1) + , _theEtaLow(-1.1) + , _thePtHigh(999.9) + , _thePtLow(0) + , _thePHigh(999.9) + , _thePLow(-999.9) + , _thePzHigh(999.9) + , _thePzLow(-999.9) + , + + _doEtaHighCut(true) + , _doEtaLowCut(true) + , _doBothEtaCut(true) + , + + _doAbsEtaHighCut(false) + , _doAbsEtaLowCut(false) + , _doBothAbsEtaCut(false) + , + + _doPtHighCut(false) + , _doPtLowCut(false) + , _doBothPtCut(false) + , + _doPHighCut(false) + , _doPLowCut(false) + , _doBothPCut(false) + , + + _doPzHighCut(false) + , _doPzLowCut(false) + , _doBothPzCut(false) +{ + if (threshold != 0) + { + _doPtLowCut = true; + _thePtLow = threshold; + } +} + +//____________________________________________________________________________.. +int HepMCParticleTrigger::process_event(PHCompositeNode* topNode) +{ + // std::cout << "HepMCParticleTrigger::process_event(PHCompositeNode *topNode) Processing Event" << std::endl; + n_evts++; + if (this->set_event_limit == true) + { // needed to keep all HepMC output at the same number of events + if (n_good >= this->goal_event_number) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + } + bool good_event{false}; + PHHepMCGenEventMap* phg = findNode::getClass(topNode, "PHHepMCGenEventMap"); + if (!phg) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + for (PHHepMCGenEventMap::ConstIter eventIter = phg->begin(); eventIter != phg->end(); ++eventIter) + { + PHHepMCGenEvent* hepev = eventIter->second; + if (!hepev) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + HepMC::GenEvent* ev = hepev->getEvent(); + if (!ev) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + good_event = isGoodEvent(ev); + if (!good_event) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + } + if (good_event) + { + n_good++; + } + return Fun4AllReturnCodes::EVENT_OK; +} +void HepMCParticleTrigger::AddParticle(int particlePid) +{ + _theParticles.push_back(particlePid); + return; +} +void HepMCParticleTrigger::AddParticles(const std::vector& particles) +{ + for (auto p : particles) + { + _theParticles.push_back(p); + } + return; +} + +void HepMCParticleTrigger::SetPtHigh(double pt) +{ + _thePtHigh = pt; + _doPtHighCut = true; + if (_doPtLowCut) + { + _doBothPtCut = true; + } + return; +} +void HepMCParticleTrigger::SetPtLow(double pt) +{ + _thePtLow = pt; + _doPtLowCut = true; + if (_doPtHighCut) + { + _doBothPtCut = true; + } + return; +} +void HepMCParticleTrigger::SetPtHighLow(double ptHigh, double ptLow) +{ + _thePtHigh = ptHigh; + _doPtHighCut = true; + _thePtLow = ptLow; + _doPtLowCut = true; + _doBothPtCut = true; + return; +} +void HepMCParticleTrigger::SetPHigh(double pt) +{ + _thePHigh = pt; + _doPHighCut = true; + if (_doPLowCut) + { + _doBothPCut = true; + } + return; +} +void HepMCParticleTrigger::SetPLow(double pt) +{ + _thePLow = pt; + _doPLowCut = true; + if (_doPHighCut) + { + _doBothPCut = true; + } + return; +} +void HepMCParticleTrigger::SetPHighLow(double ptHigh, double ptLow) +{ + _thePHigh = ptHigh; + _doPHighCut = true; + _thePLow = ptLow; + _doPLowCut = true; + _doBothPCut = true; + return; +} +void HepMCParticleTrigger::SetPzHigh(double pt) +{ + _thePzHigh = pt; + _doPzHighCut = true; + if (_doPzLowCut) + { + _doBothPzCut = true; + } + return; +} +void HepMCParticleTrigger::SetPzLow(double pt) +{ + _thePzLow = pt; + _doPzLowCut = true; + if (_doPzHighCut) + { + _doBothPzCut = true; + } + return; +} +void HepMCParticleTrigger::SetPzHighLow(double ptHigh, double ptLow) +{ + _thePzHigh = ptHigh; + _doPzHighCut = true; + _thePzLow = ptLow; + _doPzLowCut = true; + _doBothPzCut = true; + return; +} +void HepMCParticleTrigger::SetEtaHigh(double pt) +{ + _theEtaHigh = pt; + _doEtaHighCut = true; + if (_doEtaLowCut) + { + _doBothEtaCut = true; + } + return; +} +void HepMCParticleTrigger::SetEtaLow(double pt) +{ + _theEtaLow = pt; + _doEtaLowCut = true; + if (_doEtaHighCut) + { + _doBothEtaCut = true; + } + return; +} +void HepMCParticleTrigger::SetEtaHighLow(double ptHigh, double ptLow) +{ + _theEtaHigh = ptHigh; + _doEtaHighCut = true; + _theEtaLow = ptLow; + _doEtaLowCut = true; + _doBothEtaCut = true; + return; +} +void HepMCParticleTrigger::SetAbsEtaHigh(double pt) +{ + _theEtaHigh = pt; + _doAbsEtaHighCut = true; + if (_doAbsEtaLowCut) + { + _doBothAbsEtaCut = true; + } + return; +} +void HepMCParticleTrigger::SetAbsEtaLow(double pt) +{ + _theEtaLow = pt; + _doAbsEtaLowCut = true; + if (_doAbsEtaHighCut) + { + _doBothAbsEtaCut = true; + } + return; +} +void HepMCParticleTrigger::SetAbsEtaHighLow(double ptHigh, double ptLow) +{ + _theEtaHigh = ptHigh; + _doAbsEtaHighCut = true; + _theEtaLow = ptLow; + _doAbsEtaLowCut = true; + _doBothAbsEtaCut = true; + return; +} +bool HepMCParticleTrigger::isGoodEvent(HepMC::GenEvent* e1) +{ + // this is really just the call to actually evaluate and return the filter + /*if (this->threshold == 0) + { + return true; + }*/ + std::vector n_trigger_particles = getParticles(e1); + for (auto ntp : n_trigger_particles) + { + if (ntp <= 0) + { + return false; // make sure all particles have at least 1 + } + } + return true; +} + +std::vector HepMCParticleTrigger::getParticles(HepMC::GenEvent* e1) +{ + std::vector n_trigger{}; + std::map particle_types; + for (HepMC::GenEvent::particle_const_iterator iter = e1->particles_begin(); iter != e1->particles_end(); ++iter) + { + if (m_doStableParticleOnly && ((*iter)->end_vertex() || (*iter)->status() != 1)) + { + continue; + } + auto p = (*iter)->momentum(); + float px = p.px(); + float py = p.py(); + float pz = p.pz(); + float p_M = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)); + float pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); + int pid = std::abs((*iter)->pdg_id()); + double eta = p.eta(); + if ((_doEtaHighCut || _doBothEtaCut) && eta > _theEtaHigh) + { + continue; + } + if ((_doEtaLowCut || _doBothEtaCut) && eta < _theEtaLow) + { + continue; + } + if ((_doAbsEtaHighCut || _doBothAbsEtaCut) && std::abs(eta) > _theEtaHigh) + { + continue; + } + if ((_doAbsEtaLowCut || _doBothAbsEtaCut) && std::abs(eta) < _theEtaLow) + { + continue; + } + if ((_doPtHighCut || _doBothPtCut) && pt > _thePtHigh) + { + continue; + } + if ((_doPtLowCut || _doBothPtCut) && pt < _thePtLow) + { + continue; + } + if ((_doPHighCut || _doBothPCut) && p_M > _thePHigh) + { + continue; + } + if ((_doPLowCut || _doBothPCut) && p_M < _thePLow) + { + continue; + } + if ((_doPzHighCut || _doBothPzCut) && pz > _thePzHigh) + { + continue; + } + if ((_doPzLowCut || _doBothPzCut) && pz < _thePzLow) + { + continue; + } + if (particle_types.contains(pid)) + { + particle_types[pid]++; + } + else + { + particle_types[pid] = 1; + } + } + n_trigger.reserve(_theParticles.size()); + for (auto p : _theParticles) + { + n_trigger.push_back(particleAboveThreshold(particle_types, p)); // make sure we have at least one of each required particle + } + return n_trigger; +} +int HepMCParticleTrigger::particleAboveThreshold(const std::map& n_particles, int trigger_particle) +{ + // search through for the number of identified trigger particles passing cuts + auto it = n_particles.find(std::abs(trigger_particle)); + if (it != n_particles.end()) + { + return it->second; + } + return 0; +} diff --git a/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.h b/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.h new file mode 100644 index 0000000000..1b77724c0a --- /dev/null +++ b/generators/Herwig/HepMCTrigger/HepMCParticleTrigger.h @@ -0,0 +1,113 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef HEPMCPARTICLETRIGGER_H +#define HEPMCPARTICLETRIGGER_H + +#include + +#include + +#include +#include +#include +#include + +class PHCompositeNode; +namespace HepMC +{ + class GenEvent; +} + +class HepMCParticleTrigger : public SubsysReco +{ + public: + HepMCParticleTrigger(float trigger_thresh = 10., int n_incom = 1000, bool up_lim = false, const std::string& name = "HepMCParticleTrigger"); + + ~HepMCParticleTrigger() override = default; + + /** Called for each event. + This is where you do the real work. + */ + int process_event(PHCompositeNode* topNode) override; + + /// Clean up internals after each event. + + /// Called at the end of each run. + + /// Called at the end of all processing. + + /// Reset + void AddParticles(const std::vector&); //exclusively take input in the form of a pdg_ids (22 for photon, primary use case) + void AddParticle(int); + + /* void AddParents(const std::string &parents); + void AddParents(int parent); + void AddParents(std::vector parents); + void AddParentspID(std::vector parents); + */ + void SetPtHigh(double); + void SetPtLow(double); + void SetPtHighLow(double, double); + + void SetPHigh(double); + void SetPLow(double); + void SetPHighLow(double, double); + + void SetEtaHigh(double); + void SetEtaLow(double); + void SetEtaHighLow(double, double); + + void SetAbsEtaHigh(double); + void SetAbsEtaLow(double); + void SetAbsEtaHighLow(double, double); + + void SetPzHigh(double); + void SetPzLow(double); + void SetPzHighLow(double, double); + + void SetStableParticleOnly(bool b) { m_doStableParticleOnly = b; } + + private: + bool isGoodEvent(HepMC::GenEvent* e1); + std::vector getParticles(HepMC::GenEvent* e1); + int particleAboveThreshold(const std::map& n_particles, int particle); + // std::vector _theParentsi {}; + std::vector _theParticles{}; + bool m_doStableParticleOnly{true}; + float threshold{0.}; + int goal_event_number{1000}; + int n_evts{0}; + int n_good{0}; + bool set_event_limit{false}; + + float _theEtaHigh{1.1}; + float _theEtaLow{-1.1}; + float _thePtHigh{999.9}; + float _thePtLow{-999.9}; + float _thePHigh{999.9}; + float _thePLow{-999.9}; + float _thePzHigh{999.9}; + float _thePzLow{-999.9}; + + bool _doEtaHighCut{true}; + bool _doEtaLowCut{true}; + bool _doBothEtaCut{true}; + + bool _doAbsEtaHighCut{false}; + bool _doAbsEtaLowCut{false}; + bool _doBothAbsEtaCut{false}; + + bool _doPtHighCut{false}; + bool _doPtLowCut{false}; + bool _doBothPtCut{false}; + + bool _doPHighCut{false}; + bool _doPLowCut{false}; + bool _doBothPCut{false}; + + bool _doPzHighCut{false}; + bool _doPzLowCut{false}; + bool _doBothPzCut{false}; +}; + +#endif // HEPMCPARTICLETRIGGER_H diff --git a/generators/Herwig/HepMCTrigger/Makefile.am b/generators/Herwig/HepMCTrigger/Makefile.am index 17ef6c4221..9807a9a78e 100644 --- a/generators/Herwig/HepMCTrigger/Makefile.am +++ b/generators/Herwig/HepMCTrigger/Makefile.am @@ -12,20 +12,31 @@ AM_LDFLAGS = \ `fastjet-config --libs` pkginclude_HEADERS = \ - HepMCJetTrigger.h + HepMCJetTrigger.h \ + HepMCParticleTrigger.h lib_LTLIBRARIES = \ - libHepMCJetTrigger.la + libHepMCJetTrigger.la \ + libHepMCParticleTrigger.la libHepMCJetTrigger_la_SOURCES = \ HepMCJetTrigger.cc +libHepMCParticleTrigger_la_SOURCES = \ + HepMCParticleTrigger.cc + libHepMCJetTrigger_la_LIBADD = \ -lphool \ -lSubsysReco \ -lfun4all \ -lphhepmc +libHepMCParticleTrigger_la_LIBADD = \ + -lphool \ + -lSubsysReco \ + -lfun4all \ + -lphhepmc + BUILT_SOURCES = testexternals.cc noinst_PROGRAMS = \ diff --git a/generators/PHPythia8/PHPythia8.cc b/generators/PHPythia8/PHPythia8.cc index 9aa14d4684..58b1f67e5d 100644 --- a/generators/PHPythia8/PHPythia8.cc +++ b/generators/PHPythia8/PHPythia8.cc @@ -33,6 +33,20 @@ #include #include // for operator<<, endl +/** + * @brief Construct a PHPythia8 generator instance and configure HepMC conversion. + * + * Initializes the Pythia8 engine using the path from the environment variable + * `PYTHIA8`, configures a HepMC::Pythia8ToHepMC converter to store process, + * PDF, and cross-section information, and sets the default embedding ID to 1. + * The constructor preserves and restores std::cout formatting around Pythia8 + * construction to avoid altering global stream state. + * + * If `PYTHIA8` is not set, an error message is printed and the Pythia8 instance + * remains uninitialized. + * + * @param name Name forwarded to the SubsysReco base class (module instance name). + */ PHPythia8::PHPythia8(const std::string &name) : SubsysReco(name) { @@ -45,8 +59,12 @@ PHPythia8::PHPythia8(const std::string &name) std::string thePath(charPath); thePath += "/xmldoc/"; + // the pythia8 ctor messes with the formatting, so we save the cout state here + // and restore it later + std::ios old_state(nullptr); + old_state.copyfmt(std::cout); m_Pythia8.reset(new Pythia8::Pythia(thePath)); - + std::cout.copyfmt(old_state); m_Pythia8ToHepMC.reset(new HepMC::Pythia8ToHepMC()); m_Pythia8ToHepMC->set_store_proc(true); m_Pythia8ToHepMC->set_store_pdf(true); @@ -55,6 +73,18 @@ PHPythia8::PHPythia8(const std::string &name) PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1 } +/** + * @brief Initialize the Pythia8 generator, configure nodes, and seed the RNG. + * + * Performs module initialization: reads an optional configuration file and any + * queued Pythia command strings, creates the required node tree under the + * provided top-level node, sets Pythia's random seed (mapped from PHRandomSeed + * into Pythia's valid range) and prints it for reproducibility, then calls + * Pythia8::init(). + * + * @param topNode Top-level PHCompositeNode under which generator nodes are created. + * @return int Fun4All return code; returns Fun4AllReturnCodes::EVENT_OK on success. + */ int PHPythia8::Init(PHCompositeNode *topNode) { if (!m_ConfigFileName.empty()) @@ -92,8 +122,15 @@ int PHPythia8::Init(PHCompositeNode *topNode) // print out seed so we can make this is reproducible std::cout << "PHPythia8 random seed: " << seed << std::endl; + +// pythia again messes with the cout formatting + std::ios old_state(nullptr); + old_state.copyfmt(std::cout); // save current state + m_Pythia8->init(); + std::cout.copyfmt(old_state); // restore state to saved state + return Fun4AllReturnCodes::EVENT_OK; } @@ -170,6 +207,9 @@ int PHPythia8::process_event(PHCompositeNode * /*topNode*/) bool passedGen = false; bool passedTrigger = false; // int genCounter = 0; +// pythia again messes with the cout formatting in its event loop + std::ios old_state(nullptr); + old_state.copyfmt(std::cout); // save current state while (!passedTrigger) { @@ -245,6 +285,7 @@ int PHPythia8::process_event(PHCompositeNode * /*topNode*/) if (!success) { std::cout << "PHPythia8::process_event - Failed to add event to HepMC record!" << std::endl; + std::cout.copyfmt(old_state); // restore state to saved state return Fun4AllReturnCodes::ABORTRUN; } @@ -265,6 +306,8 @@ int PHPythia8::process_event(PHCompositeNode * /*topNode*/) ++m_EventCount; + std::cout.copyfmt(old_state); // restore state to saved state + // save statistics if (m_IntegralNode) { diff --git a/generators/phhepmc/PHHepMCGenEventv1.cc b/generators/phhepmc/PHHepMCGenEventv1.cc index bc8eef6f41..3aa7515365 100644 --- a/generators/phhepmc/PHHepMCGenEventv1.cc +++ b/generators/phhepmc/PHHepMCGenEventv1.cc @@ -9,6 +9,7 @@ #include #include // for cout +#include #include // for map #include #include // for swap @@ -109,6 +110,6 @@ float PHHepMCGenEventv1::get_flow_psi(unsigned int n) const return it->second; } - std::cout << "PHHepMCGenEventv1::get_flow_psi - Warning - requested reaction plane angle psi_n for n=" << n << " does not exist. Returning 0.0" << std::endl; - return 0.0F; + std::cout << "PHHepMCGenEventv1::get_flow_psi - Warning - requested reaction plane angle psi_n for n=" << n << " does not exist. Returning NAN" << std::endl; + return std::numeric_limits::quiet_NaN(); } diff --git a/offline/QA/Calorimeters/CaloValid.cc b/offline/QA/Calorimeters/CaloValid.cc index 2451672807..19c834bdab 100644 --- a/offline/QA/Calorimeters/CaloValid.cc +++ b/offline/QA/Calorimeters/CaloValid.cc @@ -335,13 +335,16 @@ int CaloValid::process_towers(PHCompositeNode* topNode) status = status >> 1U; // clang-tidy mark 1 as unsigned } - totalcemc += offlineenergy; + if (isGood) + { + totalcemc += offlineenergy; + } h_emcaltime->Fill(_timef); if (offlineenergy > emcal_hit_threshold) { h_cemc_etaphi_time->Fill(ieta, iphi, _timef); h_cemc_etaphi->Fill(ieta, iphi); - if (isGood && (scaledBits[10] || scaledBits[11])) + if (isGood && (scaledBits[10] || scaledBits[12])) { h_cemc_etaphi_wQA->Fill(ieta, iphi, offlineenergy); } @@ -407,14 +410,17 @@ int CaloValid::process_towers(PHCompositeNode* topNode) status = status >> 1U; // clang-tidy mark 1 as unsigned } - totalihcal += offlineenergy; + if (isGood) + { + totalihcal += offlineenergy; + } h_ihcaltime->Fill(_timef); if (offlineenergy > ihcal_hit_threshold) { h_ihcal_etaphi->Fill(ieta, iphi); h_ihcal_etaphi_time->Fill(ieta, iphi, _timef); - if (isGood && (scaledBits[10] || scaledBits[11])) + if (isGood && (scaledBits[10] || scaledBits[12])) { h_ihcal_etaphi_wQA->Fill(ieta, iphi, offlineenergy); } @@ -472,14 +478,17 @@ int CaloValid::process_towers(PHCompositeNode* topNode) status = status >> 1U; // clang-tidy mark 1 as unsigned } - totalohcal += offlineenergy; + if (isGood) + { + totalohcal += offlineenergy; + } h_ohcaltime->Fill(_timef); if (offlineenergy > ohcal_hit_threshold) { h_ohcal_etaphi_time->Fill(ieta, iphi, _timef); h_ohcal_etaphi->Fill(ieta, iphi); - if (isGood && (scaledBits[10] || scaledBits[11])) + if (isGood && (scaledBits[10] || scaledBits[12])) { h_ohcal_etaphi_wQA->Fill(ieta, iphi, offlineenergy); } diff --git a/offline/QA/Calorimeters/CaloValid.h b/offline/QA/Calorimeters/CaloValid.h index 2e2f9e278e..baaede0fe5 100644 --- a/offline/QA/Calorimeters/CaloValid.h +++ b/offline/QA/Calorimeters/CaloValid.h @@ -55,7 +55,7 @@ class CaloValid : public SubsysReco TriggerAnalyzer* trigAna{nullptr}; TH3* h_pi0_trigIB_mass{nullptr}; - std::vector triggerIndices{10, 28, 29, 30, 31}; // MBD NS>=1, Photon Triggers + std::vector triggerIndices{10, 12, 28, 29, 30, 31}; // MBD NS>=1, Photon Triggers TH1* h_cemc_channel_pedestal[128 * 192]{nullptr}; TH1* h_ihcal_channel_pedestal[32 * 48]{nullptr}; diff --git a/offline/QA/Tracking/Makefile.am b/offline/QA/Tracking/Makefile.am index f5d87546f1..75cbdf164b 100644 --- a/offline/QA/Tracking/Makefile.am +++ b/offline/QA/Tracking/Makefile.am @@ -17,6 +17,7 @@ pkginclude_HEADERS = \ TpcSeedsQA.h \ TpcSiliconQA.h \ SiliconSeedsQA.h \ + StateClusterResidualsQA.h \ MicromegasClusterQA.h \ CosmicTrackQA.h \ TrackFittingQA.h \ @@ -32,6 +33,7 @@ libtrackingqa_la_SOURCES = \ TpcSeedsQA.cc \ TpcSiliconQA.cc \ SiliconSeedsQA.cc \ + StateClusterResidualsQA.cc \ MicromegasClusterQA.cc \ CosmicTrackQA.cc \ TrackFittingQA.cc \ diff --git a/offline/QA/Tracking/MicromegasClusterQA.cc b/offline/QA/Tracking/MicromegasClusterQA.cc index 41efeffaa7..6feb510c54 100644 --- a/offline/QA/Tracking/MicromegasClusterQA.cc +++ b/offline/QA/Tracking/MicromegasClusterQA.cc @@ -73,6 +73,9 @@ int MicromegasClusterQA::InitRun(PHCompositeNode* topNode) << (m_calibration_filename.empty() ? "unspecified" : m_calibration_filename) << std::endl; + std::cout << "MicromegasClusterQA::InitRun - m_sample_min: " << m_sample_min << std::endl; + std::cout << "MicromegasClusterQA::InitRun - m_sample_max: " << m_sample_max << std::endl; + // read calibrations if (!m_calibration_filename.empty()) { @@ -162,6 +165,13 @@ int MicromegasClusterQA::process_event(PHCompositeNode* topNode) // find associated hits const auto hit_range = m_cluster_hit_map->getHits(ckey); + // check hit samples + // if none of the associated hits' sample is within acceptable range, skip the cluster + if( std::none_of( hit_range.first, hit_range.second, + [this]( const TrkrClusterHitAssoc::Map::value_type& pair ) + { return MicromegasDefs::getSample( pair.second ) >= m_sample_min && MicromegasDefs::getSample( pair.second ) < m_sample_max; } ) ) + { continue; } + // store cluster size and fill cluster size histogram const int cluster_size = std::distance(hit_range.first, hit_range.second); m_h_cluster_size->Fill(detid, cluster_size); diff --git a/offline/QA/Tracking/MicromegasClusterQA.h b/offline/QA/Tracking/MicromegasClusterQA.h index 338e5302de..d6de8344c4 100644 --- a/offline/QA/Tracking/MicromegasClusterQA.h +++ b/offline/QA/Tracking/MicromegasClusterQA.h @@ -51,6 +51,13 @@ class MicromegasClusterQA : public SubsysReco m_calibration_filename = value; } + /// set min sample for signal hits + void set_sample_min(uint16_t value) { m_sample_min = value; } + + /// set max sample for signal hits + void set_sample_max(uint16_t value) { m_sample_max = value; } + + private: void create_histograms(); @@ -98,6 +105,12 @@ class MicromegasClusterQA : public SubsysReco /// keep track of detector names std::vector m_detector_names; + /// min sample for signal + uint16_t m_sample_min = 0; + + /// max sample for signal + uint16_t m_sample_max = 1024; + ///@name calibration filename //@{ diff --git a/offline/QA/Tracking/StateClusterResidualsQA.cc b/offline/QA/Tracking/StateClusterResidualsQA.cc new file mode 100644 index 0000000000..1ca70dc294 --- /dev/null +++ b/offline/QA/Tracking/StateClusterResidualsQA.cc @@ -0,0 +1,245 @@ +#include "StateClusterResidualsQA.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +namespace +{ + template + inline T square (T const& t) { return t * t; } + + template + class range_adaptor + { + public: + explicit range_adaptor( + T const& begin, + T const& end) + : m_begin(begin) + , m_end(end) + { + } + T const& begin() { return m_begin; } + T const& end() { return m_end; } + + private: + T m_begin; + T m_end; + }; +} // namespace + +StateClusterResidualsQA::StateClusterResidualsQA(const std::string& name) + : SubsysReco(name) +{ +} + +int StateClusterResidualsQA::InitRun( + PHCompositeNode* top_node) +{ + createHistos(); + + // F4A will not actually ABORTRUN unless that return code is issued here + auto* track_map = findNode::getClass(top_node, m_track_map_node_name); + if (!track_map) + { + std::cout + << PHWHERE << "\n" + << "\tCould not get track map:\n" + << "\t\"" << m_track_map_node_name << "\"\n" + << "\tAborting\n" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + auto* cluster_map = findNode::getClass(top_node, m_clusterContainerName); + if (!cluster_map) + { + std::cout + << PHWHERE << "\n" + << "\tCould not get cluster map:\n" + << "\t\"" << m_clusterContainerName << "\"\n" + << "\tAborting\n" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + auto *geometry = findNode::getClass(top_node, "ActsGeometry"); + if (!geometry) + { + std::cout + << PHWHERE << "\n" + << "\tCould not get ActsGeometry:\n" + << "\t\"" << "ActsGeometry" << "\"\n" + << "\tAborting\n" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + auto* hm = QAHistManagerDef::getHistoManager(); + if (!hm) + { + std::cout + << PHWHERE << "\n" + << "\tCould not get QAHistManager\n" + << "\tAborting\n" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + for (const auto& cfg : m_pending) + { + m_histograms_x.push_back(dynamic_cast(hm->getHisto(std::string(cfg.name + "_x")))); + m_histograms_y.push_back(dynamic_cast(hm->getHisto(std::string(cfg.name + "_y")))); + m_histograms_z.push_back(dynamic_cast(hm->getHisto(std::string(cfg.name + "_z")))); + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +int StateClusterResidualsQA::process_event(PHCompositeNode* top_node) +{ + auto* track_map = findNode::getClass(top_node, m_track_map_node_name); + auto *cluster_map = findNode::getClass(top_node, m_clusterContainerName); + auto *geometry = findNode::getClass(top_node, "ActsGeometry"); + + for (auto const& [idkey, track] : *track_map) + { + if (!track) + { + continue; + } + + // count states + std::map counters = { + {TrkrDefs::mvtxId, 0}, + {TrkrDefs::inttId, 0}, + {TrkrDefs::tpcId, 0}, + {TrkrDefs::micromegasId, 0}, + }; + + for (auto const& [path_length, state] : range_adaptor(track->begin_states(), track->end_states())) + { + // There is an additional state representing the vertex at the beginning of the map, + // but getTrkrId will return 0 for its corresponding cluster + // Identify it as having path_length identically equal to 0 + if (path_length == 0) { continue; } + + auto trkr_id = static_cast(TrkrDefs::getTrkrId(state->get_cluskey())); + auto itr = counters.find(trkr_id); + if (itr == counters.end()) { continue; } + ++itr->second; + } + + float track_eta = track->get_eta(); + float track_phi = track->get_phi(); + float track_pt = track->get_pt(); + int h = 0; + for (const auto& cfg : m_pending) + { + if (cfg.charge != 0) + { + if ((cfg.charge < 0) && track->get_positive_charge()) + { + continue; + } + if ((cfg.charge > 0) && !(track->get_positive_charge())) + { + continue; + } + } + if (cfg.min_mvtx_clusters <= counters[TrkrDefs::mvtxId] && cfg.max_mvtx_clusters >= counters[TrkrDefs::mvtxId] + && cfg.min_intt_clusters <= counters[TrkrDefs::inttId] && cfg.max_intt_clusters >= counters[TrkrDefs::inttId] + && cfg.min_tpc_clusters <= counters[TrkrDefs::tpcId] && cfg.max_tpc_clusters >= counters[TrkrDefs::tpcId] + && cfg.phi_min <= track_phi && cfg.phi_max >= track_phi + && cfg.eta_min <= track_eta && cfg.eta_max >= track_eta + && cfg.pt_min <= track_pt && cfg.pt_max >= track_pt) + { + for (auto const& [path_length, state] : range_adaptor(track->begin_states(), track->end_states())) + { + if (path_length == 0) { continue; } + + auto *cluster = cluster_map->findCluster(state->get_cluskey()); + float state_x = state->get_x(); + float state_y = state->get_y(); + float state_z = state->get_z(); + Acts::Vector3 glob = geometry->getGlobalPosition(state->get_cluskey(), cluster); + float cluster_x = glob.x(); + float cluster_y = glob.y(); + float cluster_z = glob.z(); + if (cluster) + { + m_histograms_x[h]->Fill(state_x - cluster_x); + m_histograms_y[h]->Fill(state_y - cluster_y); + m_histograms_z[h]->Fill(state_z - cluster_z); + } + } + } + ++h; + } + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +void StateClusterResidualsQA::createHistos() +{ + auto *hm = QAHistManagerDef::getHistoManager(); + assert(hm); + + for (const auto& cfg : m_pending) + { + TH1F* h_new_x = new TH1F( + (cfg.name + "_x").c_str(), + ";State-Cluster X Residual [cm];Entries", + m_nBins, m_xrange.first, m_xrange.second); + h_new_x->SetMarkerColor(kBlue); + h_new_x->SetLineColor(kBlue); + hm->registerHisto(h_new_x); + TH1F* h_new_y = new TH1F( + (cfg.name + "_y").c_str(), + ";State-Cluster Y Residual [cm];Entries", + m_nBins, m_yrange.first, m_yrange.second); + h_new_y->SetMarkerColor(kBlue); + h_new_y->SetLineColor(kBlue); + hm->registerHisto(h_new_y); + TH1F* h_new_z = new TH1F( + (cfg.name + "_z").c_str(), + ";State-Cluster Z Residual [cm];Entries", + m_nBins, m_zrange.first, m_zrange.second); + h_new_z->SetMarkerColor(kBlue); + h_new_z->SetLineColor(kBlue); + hm->registerHisto(h_new_z); + } +} + +int StateClusterResidualsQA::EndRun(const int /*unused*/) +{ + auto *hm = QAHistManagerDef::getHistoManager(); + assert(hm); + + return Fun4AllReturnCodes::EVENT_OK; +} diff --git a/offline/QA/Tracking/StateClusterResidualsQA.h b/offline/QA/Tracking/StateClusterResidualsQA.h new file mode 100644 index 0000000000..739c816609 --- /dev/null +++ b/offline/QA/Tracking/StateClusterResidualsQA.h @@ -0,0 +1,127 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef STATECLUSTERRESIDUALSQA_H +#define STATECLUSTERRESIDUALSQA_H + +#include + +#include +#include +#include +#include +#include + +class PHCompositeNode; +class TH1; + +struct ResidualHistConfig +{ + std::string name = "h_StateClusterResidualsQA_"; + std::string title = ";Residual [cm];Entries"; + + int min_mvtx_clusters = 0; + int max_mvtx_clusters = 3; + int min_intt_clusters = 0; + int max_intt_clusters = 4; + int min_tpc_clusters = 0; + int max_tpc_clusters = 48; + + float phi_min = -M_PI; + float phi_max = M_PI; + float eta_min = -1.1; + float eta_max = 1.1; + + float pt_min = 0.0; + float pt_max = FLT_MAX; + + int charge = 0; +}; + +class StateClusterResidualsQA : public SubsysReco +{ + public: + StateClusterResidualsQA(const std::string& name = "StateClusterResidualsQA"); + ~StateClusterResidualsQA() override = default; + + /// sets the name of node to retrieve the track map from (default member value is "SvtxTrackMap") + void set_track_map_name(std::string const& track_map_node_name) { m_track_map_node_name = track_map_node_name; } + + StateClusterResidualsQA& addHistogram(const std::string& name) + { + ResidualHistConfig cfg; + cfg.name += name; + m_pending.push_back(cfg); + return *this; + } + StateClusterResidualsQA& setNMvtx(int min, int max) + { + m_pending.back().min_mvtx_clusters = min; + m_pending.back().max_mvtx_clusters = max; + return *this; + } + StateClusterResidualsQA& setNIntt(int min, int max) + { + m_pending.back().min_intt_clusters = min; + m_pending.back().max_intt_clusters = max; + return *this; + } + StateClusterResidualsQA& setNTpc(int min, int max) + { + m_pending.back().min_tpc_clusters = min; + m_pending.back().max_tpc_clusters = max; + return *this; + } + StateClusterResidualsQA& setPhiRange(float min, float max) + { + m_pending.back().phi_min = min; + m_pending.back().phi_max = max; + return *this; + } + StateClusterResidualsQA& setEtaRange(float min, float max) + { + m_pending.back().eta_min = min; + m_pending.back().eta_max = max; + return *this; + } + StateClusterResidualsQA& setPtRange(float min, float max) + { + m_pending.back().pt_min = min; + m_pending.back().pt_max = max; + return *this; + } + StateClusterResidualsQA& setPositiveTracks() + { + m_pending.back().charge = 1; + return *this; + } + StateClusterResidualsQA& setNegativeTracks() + { + m_pending.back().charge = -1; + return *this; + } + + void createHistos(); + + int InitRun(PHCompositeNode*) override; + + int process_event(PHCompositeNode*) override; + + int EndRun(const int runnumber) override; + + private: + std::vector m_pending; + + std::string m_track_map_node_name = "SvtxTrackMap"; + std::string m_clusterContainerName = "TRKR_CLUSTER"; + + int m_nBins = 50; + std::pair m_xrange {-0.5,0.5}; + std::pair m_yrange {-0.5,0.5}; + std::pair m_zrange {-0.5,0.5}; + + std::vector m_histograms_x{}; + std::vector m_histograms_y{}; + std::vector m_histograms_z{}; +}; + +#endif // TRACKFITTINGQA_H diff --git a/offline/framework/ffamodules/CDBInterface.cc b/offline/framework/ffamodules/CDBInterface.cc index 2b4feef203..2a0b3e1b3c 100644 --- a/offline/framework/ffamodules/CDBInterface.cc +++ b/offline/framework/ffamodules/CDBInterface.cc @@ -185,11 +185,14 @@ std::string CDBInterface::getUrl(const std::string &domain, const std::string &f std::cout << "... reply: " << return_url << std::endl; } } - auto pret = m_UrlVector.insert(make_tuple(domain_noconst, return_url, timestamp)); - if (!pret.second && Verbosity() > 1) + if (! return_url.empty()) { - std::cout << PHWHERE << "not adding again " << domain_noconst << ", url: " << return_url - << ", time stamp: " << timestamp << std::endl; + auto pret = m_UrlVector.insert(make_tuple(domain_noconst, return_url, timestamp)); + if (!pret.second && Verbosity() > 1) + { + std::cout << PHWHERE << "not adding again " << domain_noconst << ", url: " << return_url + << ", time stamp: " << timestamp << std::endl; + } } return return_url; } diff --git a/offline/framework/ffamodules/CDBInterface.h b/offline/framework/ffamodules/CDBInterface.h index c88fae9dee..4132b35a34 100644 --- a/offline/framework/ffamodules/CDBInterface.h +++ b/offline/framework/ffamodules/CDBInterface.h @@ -10,7 +10,6 @@ #include #include // for tuple -class PHCompositeNode; class SphenixClient; class CDBInterface : public SubsysReco diff --git a/offline/framework/ffamodules/FlagHandler.h b/offline/framework/ffamodules/FlagHandler.h index f5cc06bcdb..fd7e40cca3 100644 --- a/offline/framework/ffamodules/FlagHandler.h +++ b/offline/framework/ffamodules/FlagHandler.h @@ -7,14 +7,12 @@ #include -class PHCompositeNode; - class FlagHandler : public SubsysReco { public: FlagHandler(const std::string &name = "FlagHandler"); - ~FlagHandler() override {} + ~FlagHandler() override = default; /** Create the Flag Node if it does not exist, if it exists, read back flags and copy them into recoConsts diff --git a/offline/framework/ffamodules/HeadReco.cc b/offline/framework/ffamodules/HeadReco.cc index 75e879fdff..17214758c8 100644 --- a/offline/framework/ffamodules/HeadReco.cc +++ b/offline/framework/ffamodules/HeadReco.cc @@ -84,10 +84,13 @@ int HeadReco::process_event(PHCompositeNode *topNode) { evtheader->set_ImpactParameter(hi->impact_parameter()); evtheader->set_EventPlaneAngle(hi->event_plane_angle()); - for (unsigned int n = 1; n <= 6; ++n) - { - evtheader->set_FlowPsiN(n, genevt->get_flow_psi(n)); - } + if (! genevt->get_flow_psi_map().empty()) + { + for (unsigned int n = 1; n <= 6; ++n) + { + evtheader->set_FlowPsiN(n, genevt->get_flow_psi(n)); + } + } evtheader->set_eccentricity(hi->eccentricity()); evtheader->set_ncoll(hi->Ncoll()); evtheader->set_npart(hi->Npart_targ() + hi->Npart_proj()); diff --git a/offline/framework/ffamodules/HeadReco.h b/offline/framework/ffamodules/HeadReco.h index ca2f6cf34e..53fbc3c6c3 100644 --- a/offline/framework/ffamodules/HeadReco.h +++ b/offline/framework/ffamodules/HeadReco.h @@ -7,13 +7,11 @@ #include // for string -class PHCompositeNode; - class HeadReco : public SubsysReco { public: HeadReco(const std::string &name = "HeadReco"); - ~HeadReco() override {} + ~HeadReco() override = default; int Init(PHCompositeNode *topNode) override; int InitRun(PHCompositeNode *topNode) override; int process_event(PHCompositeNode *topNode) override; diff --git a/offline/framework/ffamodules/SyncReco.h b/offline/framework/ffamodules/SyncReco.h index 0a490302a0..73d0954593 100644 --- a/offline/framework/ffamodules/SyncReco.h +++ b/offline/framework/ffamodules/SyncReco.h @@ -5,13 +5,11 @@ #include -class PHCompositeNode; - class SyncReco : public SubsysReco { public: SyncReco(const std::string &name = "SYNC"); - ~SyncReco() override {} + ~SyncReco() override = default; int Init(PHCompositeNode *topNode) override; int InitRun(PHCompositeNode *topNode) override; @@ -24,7 +22,7 @@ class SyncReco : public SubsysReco // just if we need to override the segment for e.g. embedding // where we want to reuse hijing files which normally set // the segment number - int forced_segment = -1; + int forced_segment {-1}; }; #endif /* FFAMODULES_SYNCRECO_H */ diff --git a/offline/framework/ffamodules/Timing.cc b/offline/framework/ffamodules/Timing.cc index f662ba1b33..16210e4f73 100644 --- a/offline/framework/ffamodules/Timing.cc +++ b/offline/framework/ffamodules/Timing.cc @@ -3,8 +3,6 @@ #include #include // for SubsysReco -#include - #include Timing::Timing(const std::string &name) diff --git a/offline/framework/ffamodules/Timing.h b/offline/framework/ffamodules/Timing.h index 01b95638e0..a9a5ccd39d 100644 --- a/offline/framework/ffamodules/Timing.h +++ b/offline/framework/ffamodules/Timing.h @@ -8,15 +8,13 @@ #include // for string #include -class PHCompositeNode; - class Timing : public SubsysReco { public: Timing(const std::string &name = "Timing"); - ~Timing() override {} - int InitRun(PHCompositeNode *topNode) override; - int process_event(PHCompositeNode *topNode) override; + ~Timing() override = default; + int InitRun(PHCompositeNode * /*topNode*/) override; + int process_event(PHCompositeNode * /*topNode*/) override; void SetCallCounter(unsigned int i) { calls = i; } private: diff --git a/offline/framework/frog/CreateFileList.pl b/offline/framework/frog/CreateFileList.pl index cc59083d8c..19b552b764 100755 --- a/offline/framework/frog/CreateFileList.pl +++ b/offline/framework/frog/CreateFileList.pl @@ -74,7 +74,8 @@ "33" => "JS pythia8 Jet ptmin = 15GeV", "34" => "JS pythia8 Jet ptmin = 50GeV", "35" => "JS pythia8 Jet ptmin = 70GeV", - "36" => "JS pythia8 Jet ptmin = 5GeV" + "36" => "JS pythia8 Jet ptmin = 5GeV", + "37" => "hijing O+O (0-15fm)" ); my %pileupdesc = ( @@ -926,6 +927,20 @@ $pileupstring = $pp_pileupstring; &commonfiletypes(); } + elsif ($prodtype == 37) + { + if (defined $nopileup) + { + $filenamestring = sprintf("sHijing_OO_0_15fm"); + } + else + { + $filenamestring = sprintf("sHijing_OO_0_15fm%s",$AuAu_pileupstring); + } + $notlike{$filenamestring} = ["pythia8" ,"single", "special"]; + $pileupstring = $AuAu_pileupstring; + &commonfiletypes(); + } else { diff --git a/offline/framework/fun4all/Fun4AllServer.cc b/offline/framework/fun4all/Fun4AllServer.cc index 9d8779204e..9dc68afbb0 100644 --- a/offline/framework/fun4all/Fun4AllServer.cc +++ b/offline/framework/fun4all/Fun4AllServer.cc @@ -128,9 +128,9 @@ void Fun4AllServer::InitAll() { gSystem->IgnoreSignal((ESignals) i); } + m_saved_cout_state.copyfmt(std::cout); // save current state Fun4AllMonitoring::instance()->Snapshot("StartUp"); - std::string histomanagername; - histomanagername = Name() + "HISTOS"; + std::string histomanagername = Name() + "HISTOS"; ServerHistoManager = new Fun4AllHistoManager(histomanagername); registerHistoManager(ServerHistoManager); double uplim = NFRAMEWORKBINS - 0.5; @@ -245,6 +245,7 @@ int Fun4AllServer::registerSubsystem(SubsysReco *subsystem, const std::string &t << subsystem->Name() << std::endl; exit(1); } + std::cout.copyfmt(m_saved_cout_state); // restore cout to default formatting gROOT->cd(currdir.c_str()); if (iret) { @@ -576,6 +577,7 @@ int Fun4AllServer::process_event() ffamemtracker->Snapshot("Fun4AllServerProcessEvent"); #endif int retcode = Subsystem.first->process_event(Subsystem.second); + std::cout.copyfmt(m_saved_cout_state); // restore cout to default formatting #ifdef FFAMEMTRACKER ffamemtracker->Snapshot("Fun4AllServerProcessEvent"); #endif @@ -899,6 +901,7 @@ int Fun4AllServer::BeginRun(const int runno) for (iter = Subsystems.begin(); iter != Subsystems.end(); ++iter) { iret = BeginRunSubsystem(*iter); + std::cout.copyfmt(m_saved_cout_state); // restore cout to default formatting } for (; !NewSubsystems.empty(); NewSubsystems.pop_front()) { @@ -1092,6 +1095,7 @@ int Fun4AllServer::EndRun(const int runno) << (*iter).first->Name() << std::endl; exit(1); } + std::cout.copyfmt(m_saved_cout_state); // restore cout to default formatting } gROOT->cd(currdir.c_str()); @@ -1144,6 +1148,7 @@ int Fun4AllServer::End() << (*iter).first->Name() << std::endl; exit(1); } + std::cout.copyfmt(m_saved_cout_state); // restore cout to default formatting } gROOT->cd(currdir.c_str()); PHNodeIterator nodeiter(TopNode); diff --git a/offline/framework/fun4all/Fun4AllServer.h b/offline/framework/fun4all/Fun4AllServer.h index 0c7ea72c14..9f360b9502 100644 --- a/offline/framework/fun4all/Fun4AllServer.h +++ b/offline/framework/fun4all/Fun4AllServer.h @@ -143,7 +143,8 @@ class Fun4AllServer : public Fun4AllBase int eventnumber{0}; int eventcounter{0}; int keep_db_connected{0}; - + + std::ios m_saved_cout_state{nullptr}; std::vector ComplaintList; std::vector ResetNodeList {"DST"}; std::vector> Subsystems; diff --git a/offline/framework/fun4all/InputFileHandler.cc b/offline/framework/fun4all/InputFileHandler.cc index 8939b67568..0fdc858db8 100644 --- a/offline/framework/fun4all/InputFileHandler.cc +++ b/offline/framework/fun4all/InputFileHandler.cc @@ -3,10 +3,13 @@ #include +#include + #include #include #include #include +#include int InputFileHandler::AddFile(const std::string &filename) { @@ -89,6 +92,19 @@ int InputFileHandler::OpenNextFile() { std::cout << PHWHERE << " opening next file: " << *iter << std::endl; } + if (!GetOpeningScript().empty()) + { + std::vector stringvec; + stringvec.push_back(*iter); + if (!m_FileName.empty()) + { + stringvec.push_back(m_FileName); + } + if (RunBeforeOpening(stringvec)) + { + std::cout << PHWHERE << " RunBeforeOpening() failed" << std::endl; + } + } if (fileopen(*iter)) { std::cout << PHWHERE << " could not open file: " << *iter << std::endl; @@ -145,3 +161,40 @@ int InputFileHandler::fileopen(const std::string &fname) std::cout << "InputFileHandler::fileopen opening " << fname << std::endl; return 0; } + +int InputFileHandler::RunBeforeOpening(const std::vector &stringvec) +{ + if (m_RunBeforeOpeningScript.empty()) + { + return 0; + } + if (!std::filesystem::exists(m_RunBeforeOpeningScript)) + { + std::cout << PHWHERE << " script " << m_RunBeforeOpeningScript << " not found" + << std::endl; + return -1; + } + if (!((std::filesystem::status(m_RunBeforeOpeningScript).permissions() & std::filesystem::perms::owner_exec) == std::filesystem::perms::owner_exec)) + { + std::cout << PHWHERE << "RunBeforeOpeningScript script " + << m_RunBeforeOpeningScript << " is not owner executable" << std::endl; + return -1; + } + std::string fullcmd = m_RunBeforeOpeningScript + " " + m_OpeningArgs; + for (const auto& iter : stringvec) + { + fullcmd += " " + iter; + } + + if (m_Verbosity > 1) + { + std::cout << PHWHERE << " running " << fullcmd << std::endl; + } + unsigned int iret = gSystem->Exec(fullcmd.c_str()); + + if (iret) + { + iret = iret >> 8U; + } + return static_cast (iret); +} diff --git a/offline/framework/fun4all/InputFileHandler.h b/offline/framework/fun4all/InputFileHandler.h index 823ae33b38..49df4aa379 100644 --- a/offline/framework/fun4all/InputFileHandler.h +++ b/offline/framework/fun4all/InputFileHandler.h @@ -4,13 +4,14 @@ #include #include #include +#include class InputFileHandler { public: InputFileHandler() = default; virtual ~InputFileHandler() = default; - virtual int fileopen(const std::string & /*filename*/);// { return 0; } + virtual int fileopen(const std::string & /*filename*/); // { return 0; } virtual int fileclose() { return -1; } virtual int ResetFileList(); @@ -32,12 +33,19 @@ class InputFileHandler std::pair::const_iterator, std::list::const_iterator> FileOpenListBeginEnd() { return std::make_pair(m_FileListOpened.begin(), m_FileListOpened.end()); } const std::list &GetFileList() const { return m_FileListCopy; } const std::list &GetFileOpenedList() const { return m_FileListOpened; } + void SetOpeningScript(const std::string &script) { m_RunBeforeOpeningScript = script; } + const std::string &GetOpeningScript() const { return m_RunBeforeOpeningScript; } + void SetOpeningScriptArgs(const std::string &args) { m_OpeningArgs = args; } + const std::string &GetOpeningScriptArgs() const { return m_OpeningArgs; } + int RunBeforeOpening(const std::vector &stringvec); private: int m_IsOpen{0}; int m_Repeat{0}; uint64_t m_Verbosity{0}; std::string m_FileName; + std::string m_RunBeforeOpeningScript; + std::string m_OpeningArgs; std::list m_FileList; std::list m_FileListCopy; std::list m_FileListOpened; // all files which were opened during running diff --git a/offline/framework/fun4allraw/Fun4AllStreamingInputManager.cc b/offline/framework/fun4allraw/Fun4AllStreamingInputManager.cc index 78b3e3a365..0bea81197b 100644 --- a/offline/framework/fun4allraw/Fun4AllStreamingInputManager.cc +++ b/offline/framework/fun4allraw/Fun4AllStreamingInputManager.cc @@ -776,9 +776,15 @@ int Fun4AllStreamingInputManager::FillIntt() { h_taggedAllFee_intt->Fill(refbcobitshift); } - while (m_InttRawHitMap.begin()->first <= select_crossings - m_intt_negative_bco) + + for (auto& [bco, hitinfo] : m_InttRawHitMap) { - for (auto *intthititer : m_InttRawHitMap.begin()->second.InttRawHitVector) + if (bco > select_crossings) + { + break; + } + + for (auto *intthititer : hitinfo.InttRawHitVector) { if (Verbosity() > 1) { @@ -788,25 +794,9 @@ int Fun4AllStreamingInputManager::FillIntt() } inttcont->AddHit(intthititer); } - for (auto *iter : m_InttInputVector) - { - iter->CleanupUsedPackets(m_InttRawHitMap.begin()->first); - if (m_intt_negative_bco < 2) // triggered mode - { - iter->clearPacketBClkStackMap(m_InttRawHitMap.begin()->first); - iter->clearFeeGTML1BCOMap(m_InttRawHitMap.begin()->first); - } - } - m_InttRawHitMap.begin()->second.InttRawHitVector.clear(); - m_InttRawHitMap.erase(m_InttRawHitMap.begin()); - if (m_InttRawHitMap.empty()) - { - break; - } - } + } return 0; } - int Fun4AllStreamingInputManager::FillMvtx() { int iret = FillMvtxPool(); @@ -846,7 +836,7 @@ int Fun4AllStreamingInputManager::FillMvtx() } select_crossings += m_RefBCO; - uint64_t ref_bco_minus_range = m_RefBCO < m_mvtx_bco_range ? 0 : m_RefBCO - m_mvtx_bco_range; + uint64_t ref_bco_minus_range = m_RefBCO < m_mvtx_negative_bco ? 0 : m_RefBCO - m_mvtx_negative_bco; if (Verbosity() > 2) { std::cout << "select MVTX crossings" @@ -981,90 +971,42 @@ int Fun4AllStreamingInputManager::FillMvtx() } taggedPacketsFEEs.clear(); - if (m_mvtx_is_triggered) + uint64_t lower_limit = m_mvtx_is_triggered ? select_crossings : select_crossings - m_mvtx_bco_range - m_mvtx_negative_bco; + uint64_t upper_limit = m_mvtx_is_triggered ? select_crossings + m_mvtx_bco_range : select_crossings; + + for (auto& [bco, hitinfo] : m_MvtxRawHitMap) { - while (select_crossings <= m_MvtxRawHitMap.begin()->first && m_MvtxRawHitMap.begin()->first <= select_crossings + m_mvtx_bco_range) // triggered + if (bco < lower_limit) { - if (Verbosity() > 2) - { - std::cout << "Adding 0x" << std::hex << m_MvtxRawHitMap.begin()->first - << " ref: 0x" << select_crossings << std::dec << std::endl; - } - for (auto *mvtxFeeIdInfo : m_MvtxRawHitMap.begin()->second.MvtxFeeIdInfoVector) - { - if (Verbosity() > 1) - { - mvtxFeeIdInfo->identify(); - } - mvtxEvtHeader->AddFeeIdInfo(mvtxFeeIdInfo); - delete mvtxFeeIdInfo; - } - m_MvtxRawHitMap.begin()->second.MvtxFeeIdInfoVector.clear(); - mvtxEvtHeader->AddL1Trg(m_MvtxRawHitMap.begin()->second.MvtxL1TrgBco); + continue; + } + if (bco > upper_limit) + { + break; + } - for (auto *mvtxhititer : m_MvtxRawHitMap.begin()->second.MvtxRawHitVector) - { - if (Verbosity() > 1) - { - mvtxhititer->identify(); - } - mvtxcont->AddHit(mvtxhititer); - } - for (auto *iter : m_MvtxInputVector) - { - iter->CleanupUsedPackets(m_MvtxRawHitMap.begin()->first); - } - m_MvtxRawHitMap.begin()->second.MvtxRawHitVector.clear(); - m_MvtxRawHitMap.begin()->second.MvtxL1TrgBco.clear(); - m_MvtxRawHitMap.erase(m_MvtxRawHitMap.begin()); - // m_MvtxRawHitMap.empty() need to be checked here since we do not call FillPoolMvtx() - if (m_MvtxRawHitMap.empty()) - { - break; - } + if (Verbosity() > 2) + { + std::cout << "Adding 0x" << std::hex << bco + << " ref: 0x" << select_crossings << std::dec << std::endl; } - } - else - { - while (select_crossings - m_mvtx_bco_range - m_mvtx_negative_bco <= m_MvtxRawHitMap.begin()->first && m_MvtxRawHitMap.begin()->first <= select_crossings) // streamed + for (auto *mvtxFeeIdInfo : hitinfo.MvtxFeeIdInfoVector) { - if (Verbosity() > 2) - { - std::cout << "Adding 0x" << std::hex << m_MvtxRawHitMap.begin()->first - << " ref: 0x" << select_crossings << std::dec << std::endl; - } - for (auto *mvtxFeeIdInfo : m_MvtxRawHitMap.begin()->second.MvtxFeeIdInfoVector) + if (Verbosity() > 1) { - if (Verbosity() > 1) - { - mvtxFeeIdInfo->identify(); - } - mvtxEvtHeader->AddFeeIdInfo(mvtxFeeIdInfo); - delete mvtxFeeIdInfo; + mvtxFeeIdInfo->identify(); } - m_MvtxRawHitMap.begin()->second.MvtxFeeIdInfoVector.clear(); - mvtxEvtHeader->AddL1Trg(m_MvtxRawHitMap.begin()->second.MvtxL1TrgBco); + mvtxEvtHeader->AddFeeIdInfo(mvtxFeeIdInfo); + } + mvtxEvtHeader->AddL1Trg(hitinfo.MvtxL1TrgBco); - for (auto *mvtxhititer : m_MvtxRawHitMap.begin()->second.MvtxRawHitVector) - { - if (Verbosity() > 1) - { - mvtxhititer->identify(); - } - mvtxcont->AddHit(mvtxhititer); - } - for (auto *iter : m_MvtxInputVector) - { - iter->CleanupUsedPackets(m_MvtxRawHitMap.begin()->first); - } - m_MvtxRawHitMap.begin()->second.MvtxRawHitVector.clear(); - m_MvtxRawHitMap.begin()->second.MvtxL1TrgBco.clear(); - m_MvtxRawHitMap.erase(m_MvtxRawHitMap.begin()); - // m_MvtxRawHitMap.empty() need to be checked here since we do not call FillPoolMvtx() - if (m_MvtxRawHitMap.empty()) + for (auto *mvtxhititer : hitinfo.MvtxRawHitVector) + { + if (Verbosity() > 1) { - break; + mvtxhititer->identify(); } + mvtxcont->AddHit(mvtxhititer); } } @@ -1422,7 +1364,7 @@ int Fun4AllStreamingInputManager::FillMicromegasPool() int Fun4AllStreamingInputManager::FillMvtxPool() { - uint64_t ref_bco_minus_range = m_RefBCO < m_mvtx_bco_range ? m_mvtx_bco_range : m_RefBCO - m_mvtx_bco_range; + uint64_t ref_bco_minus_range = m_RefBCO < m_mvtx_negative_bco ? m_mvtx_negative_bco : m_RefBCO - m_mvtx_negative_bco; for (auto *iter : m_MvtxInputVector) { if (Verbosity() > 3) diff --git a/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.cc b/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.cc index 7500e3badf..89c0caa9f3 100644 --- a/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.cc +++ b/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.cc @@ -607,6 +607,7 @@ void SingleMicromegasPoolInput_v2::createQAHistos() m_evaluation_file.reset(new TFile(m_evaluation_filename.c_str(), "RECREATE")); m_evaluation_tree = new TTree("T", "T"); m_evaluation_tree->Branch("is_heartbeat", &m_waveform.is_heartbeat); + m_evaluation_tree->Branch("matched", &m_waveform.matched); m_evaluation_tree->Branch("packet_id", &m_waveform.packet_id); m_evaluation_tree->Branch("fee_id", &m_waveform.fee_id); m_evaluation_tree->Branch("channel", &m_waveform.channel); @@ -912,6 +913,13 @@ void SingleMicromegasPoolInput_v2::process_fee_data(int packet_id, unsigned int // try get gtm bco matching fee const auto& fee_bco = payload.bx_timestamp; + if( m_do_evaluation ) + { + m_waveform.is_heartbeat = is_heartbeat; + m_waveform.fee_id = fee_id; + m_waveform.channel = payload.channel; + m_waveform.fee_bco = fee_bco; + } // find matching gtm bco uint64_t gtm_bco = 0; @@ -920,9 +928,20 @@ void SingleMicromegasPoolInput_v2::process_fee_data(int packet_id, unsigned int { // assign gtm bco gtm_bco = result.value(); - } - else - { + if( m_do_evaluation ) + { + m_waveform.matched = true; + m_waveform.gtm_bco_matched = gtm_bco; + { + const auto predicted = bco_matching_information.get_predicted_fee_bco(gtm_bco);; + if( predicted ) + { + m_waveform.fee_bco_predicted_matched = predicted.value(); + } + } + m_evaluation_tree->Fill(); + } + } else { // increment counter and histogram ++m_waveform_counters[packet_id].dropped_bco; ++m_fee_waveform_counters[fee_id].dropped_bco; @@ -935,28 +954,16 @@ void SingleMicromegasPoolInput_v2::process_fee_data(int packet_id, unsigned int ++m_fee_heartbeat_counters[fee_id].dropped_bco; } - // skip the waverform - continue; - } - - if (m_do_evaluation) - { - m_waveform.is_heartbeat = (payload.type == HEARTBEAT_T); - m_waveform.fee_id = fee_id; - m_waveform.channel = payload.channel; - m_waveform.fee_bco = fee_bco; - - m_waveform.gtm_bco_matched = gtm_bco; + if( m_do_evaluation ) { - const auto predicted = bco_matching_information.get_predicted_fee_bco(gtm_bco); - ; - if (predicted) - { - m_waveform.fee_bco_predicted_matched = predicted.value(); - } + m_waveform.matched = false; + m_waveform.gtm_bco_matched = 0; + m_waveform.fee_bco_predicted_matched = 0; + m_evaluation_tree->Fill(); } - m_evaluation_tree->Fill(); + // skip the waverform + continue; } // ignore heartbeat waveforms diff --git a/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.h b/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.h index d8c3c1da62..3746b9eb41 100644 --- a/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.h +++ b/offline/framework/fun4allraw/SingleMicromegasPoolInput_v2.h @@ -230,6 +230,9 @@ class SingleMicromegasPoolInput_v2 : public SingleStreamingInput /// true if measurement is hearbeat bool is_heartbeat = false; + /// true if matched + bool matched = false; + /// ll1 bco uint64_t gtm_bco_first {0}; diff --git a/offline/framework/fun4allraw/SingleMvtxPoolInput.cc b/offline/framework/fun4allraw/SingleMvtxPoolInput.cc index de98a9936e..4ce498ee9e 100644 --- a/offline/framework/fun4allraw/SingleMvtxPoolInput.cc +++ b/offline/framework/fun4allraw/SingleMvtxPoolInput.cc @@ -1,7 +1,7 @@ #include "SingleMvtxPoolInput.h" -#include "MvtxRawDefs.h" #include "Fun4AllStreamingInputManager.h" +#include "MvtxRawDefs.h" #include "mvtx_pool.h" #include @@ -29,7 +29,8 @@ #include SingleMvtxPoolInput::SingleMvtxPoolInput(const std::string &name) - : SingleStreamingInput(name), plist(new Packet *[2]) + : SingleStreamingInput(name) + , plist(new Packet *[2]) { m_rawHitContainerName = "MVTXRAWHIT"; @@ -161,7 +162,7 @@ void SingleMvtxPoolInput::FillPool(const uint64_t minBCO) m_BclkStack.insert(strb_bco); m_FEEBclkMap[feeId] = strb_bco; - if (strb_bco < minBCO - m_NegativeBco) + if (strb_bco < minBCO) { continue; } @@ -206,7 +207,7 @@ void SingleMvtxPoolInput::FillPool(const uint64_t minBCO) auto it = m_BclkStack.lower_bound(lv1Bco); // auto const strb_it = (it == m_BclkStack.begin()) ? (*it == lv1Bco ? it : m_BclkStack.cend()) : --it; // this is equivalent but human readable for the above: - auto strb_it = m_BclkStack.cend(); + auto strb_it = m_BclkStack.cend(); if (it == m_BclkStack.begin()) { @@ -462,7 +463,7 @@ void SingleMvtxPoolInput::ConfigureStreamingInputManager() else if (m_strobeWidth > 9 && m_strobeWidth < 11) { m_BcoRange = 500; - m_NegativeBco = 500; + m_NegativeBco = 120; } else if (m_strobeWidth < 1) // triggered mode { diff --git a/offline/framework/fun4allraw/SingleTriggeredInput.cc b/offline/framework/fun4allraw/SingleTriggeredInput.cc index 3bb6b30a2b..80e9f49a64 100644 --- a/offline/framework/fun4allraw/SingleTriggeredInput.cc +++ b/offline/framework/fun4allraw/SingleTriggeredInput.cc @@ -352,6 +352,12 @@ int SingleTriggeredInput::FillEventVector() m_bclkarray_map[pid][0] = tmp; m_bclkdiffarray_map[pid].fill(std::numeric_limits::max()); + static bool firstclockarray=true; + if(firstclockarray){ + std::cout << "first clock call pid " << pid << " m_bclkarray_map[pid][0] : " << m_bclkarray_map[pid][0] << std::endl; + firstclockarray=false; + } + if ( representative_pid == -1 ) { representative_pid = pid; @@ -368,12 +374,17 @@ int SingleTriggeredInput::FillEventVector() while (i < pooldepth) { Event* evt{nullptr}; + bool skiptrace = false; if (this != Gl1Input()) { auto* gl1 = dynamic_cast(Gl1Input()); if (gl1) { int nskip = gl1->GetGl1SkipArray()[i]; + if(nskip >0) + { + skiptrace = true; + } while (nskip > 0) { @@ -391,6 +402,7 @@ int SingleTriggeredInput::FillEventVector() if (skip_evt->getEvtType() != DATAEVENT) { + delete skip_evt; continue; } @@ -412,14 +424,91 @@ int SingleTriggeredInput::FillEventVector() { if (Verbosity() > 0) { - std::cout << Name() << ": Early stop of SEB skip after " << (gl1->GetGl1SkipArray()[i] - nskip) << " from intial " << gl1->GetGl1SkipArray()[i] << " events." << std::endl; + std::cout << Name() << ": Early stop in pool " << i << " of SEB skip after " << (gl1->GetGl1SkipArray()[i] - nskip) << " from intial " << gl1->GetGl1SkipArray()[i] << " events. gl1diff vs sebdiff : " << gl1_diff << " vs " << seb_diff << std::endl; } evt = skip_evt; + skiptrace = false; break; } delete skip_evt; nskip--; } + + if(skiptrace) + { + evt = GetEventIterator()->getNextEvent(); + while (!evt) + { + fileclose(); + if (OpenNextFile() == InputFileHandlerReturnCodes::FAILURE) + { + FilesDone(1); + return -1; + } + evt = GetEventIterator()->getNextEvent(); + } + if (evt->getEvtType() != DATAEVENT) + { + if (Verbosity() > 0) + { + std::cout << Name() << " dropping non data event: " << evt->getEvtSequence() << std::endl; + } + delete evt; + continue; + } + + Packet* pkt = evt->getPacket(representative_pid); + if (!pkt) + { + std::cout << "representative packet invalid inside skiptrace.. continuing.." << std::endl; + continue; + } + FillPacketClock(evt, pkt, i); + uint64_t seb_diff = m_bclkdiffarray_map[representative_pid][i]; + int gl1pid = Gl1Input()->m_bclkdiffarray_map.begin()->first; + uint64_t gl1_diff = gl1->m_bclkdiffarray_map[gl1pid][i]; + + bool clockconsistency=true; + if (seb_diff != gl1_diff) + { + clockconsistency=false; + int clockconstcount = 0; + while(!clockconsistency && clockconstcount<5) + { + std::cout << Name() << ": Still inconsistent clock diff after Gl1 drop. gl1diff vs sebdiff : " << gl1_diff << " vs " << seb_diff << std::endl; + delete pkt; + delete evt; + evt = GetEventIterator()->getNextEvent(); + while (!evt) + { + fileclose(); + if (OpenNextFile() == InputFileHandlerReturnCodes::FAILURE) + { + FilesDone(1); + return -1; + } + evt = GetEventIterator()->getNextEvent(); + } + pkt = evt->getPacket(representative_pid); + if (!pkt) + { + std::cout << "representative packet invalid inside skiptrace.. continuing.." << std::endl; + continue; + } + + FillPacketClock(evt, pkt, i); + uint64_t seb_diff_next = m_bclkdiffarray_map[representative_pid][i]; + uint64_t gl1_diff_next = gl1->m_bclkdiffarray_map[gl1pid][i]; + std::cout << "seb_diff_next : " << seb_diff_next << " , gl1_diff_next : " << gl1_diff_next << std::endl; + if(seb_diff_next == gl1_diff_next) + { + clockconsistency=true; + std::cout << Name() << " : recovered by additional skip in skiptrace" << std::endl; + } + clockconstcount++; + } + } + } } } @@ -485,14 +574,16 @@ int SingleTriggeredInput::FillEventVector() } FillPacketClock(thisevt, pkt, i); m_PacketEventDeque[pid].push_back(thisevt); + delete pkt; - + if (representative_pid == -1 && m_PacketShiftOffset[pid] == 0) { representative_pid = pid; } } i++; + eventcounter++; } size_t minSize = pooldepth; @@ -1116,16 +1207,15 @@ int SingleTriggeredInput::ReadEvent() [](const std::pair& p) { return p.second == 0; }); std::set events_to_delete; - for (auto& [pid, dq] : m_PacketEventDeque) { if(m_PacketAlignmentProblem[pid]) { continue; } + Event* evt = dq.front(); Packet* packet = evt->getPacket(pid); - int packet_id = packet->getIdentifier(); if (packet_id != pid) { @@ -1137,7 +1227,6 @@ int SingleTriggeredInput::ReadEvent() CaloPacket *newhit = findNode::getClass(m_topNode, packet_id); newhit->Reset(); - if (m_DitchPackets.contains(packet_id) && m_DitchPackets[packet_id].contains(0)) { newhit->setStatus(OfflinePacket::PACKET_DROPPED); diff --git a/offline/framework/fun4allraw/SingleTriggeredInput.h b/offline/framework/fun4allraw/SingleTriggeredInput.h index 088be73bef..d7d825b02a 100644 --- a/offline/framework/fun4allraw/SingleTriggeredInput.h +++ b/offline/framework/fun4allraw/SingleTriggeredInput.h @@ -34,9 +34,6 @@ class SingleTriggeredInput : public Fun4AllBase, public InputFileHandler virtual void FillPool(); virtual void RunNumber(const int runno) { m_RunNumber = runno; } virtual int RunNumber() const { return m_RunNumber; } - virtual void EventNumber(const int i) { m_EventNumber = i; } - virtual int EventNumber() const { return m_EventNumber; } - virtual int EventsInThisFile() const { return m_EventsThisFile; } virtual int fileopen(const std::string &filename) override; virtual int fileclose() override; virtual int AllDone() const { return m_AllDone; } @@ -45,6 +42,8 @@ class SingleTriggeredInput : public Fun4AllBase, public InputFileHandler virtual void FilesDone(const int i) { m_FilesDone = i; } virtual void EventAlignmentProblem(const int i) { m_EventAlignmentProblem = i; } virtual int EventAlignmentProblem() const { return m_EventAlignmentProblem; } + virtual void EventNumber(const int i) { m_EventNumber = i; } + virtual int EventNumber() const { return m_EventNumber; } virtual void CreateDSTNodes(Event *evt); // these ones are used directly by the derived classes, maybe later // move to cleaner accessors @@ -94,7 +93,6 @@ class SingleTriggeredInput : public Fun4AllBase, public InputFileHandler int m_AllDone{0}; uint64_t m_Event{0}; int m_EventNumber{0}; - int m_EventsThisFile{0}; int m_EventAlignmentProblem{0}; int m_FilesDone{0}; int m_LastEvent{std::numeric_limits::max()}; @@ -114,6 +112,7 @@ class SingleTriggeredInput : public Fun4AllBase, public InputFileHandler std::map m_PacketAlignmentProblem; std::map m_PrevPoolLastDiffBad; std::map m_PreviousValidBCOMap; + long long eventcounter{0}; }; #endif diff --git a/offline/packages/CaloEmbedding/CombineTowerInfo.cc b/offline/packages/CaloEmbedding/CombineTowerInfo.cc new file mode 100644 index 0000000000..3e5150d194 --- /dev/null +++ b/offline/packages/CaloEmbedding/CombineTowerInfo.cc @@ -0,0 +1,99 @@ +#include "CombineTowerInfo.h" + +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +//____________________________________________________________________________ +CombineTowerInfo::CombineTowerInfo(const std::string& name) + : SubsysReco(name) +{ +} + +//____________________________________________________________________________ +int CombineTowerInfo::InitRun(PHCompositeNode* topNode) +{ + if (m_inputNodeA.empty() || m_inputNodeB.empty() || m_outputNode.empty()) + { + throw std::runtime_error("CombineTowerInfo: input/output node names not set"); + } + + CreateNodes(topNode); + return Fun4AllReturnCodes::EVENT_OK; +} + +//____________________________________________________________________________ +void CombineTowerInfo::CreateNodes(PHCompositeNode* topNode) +{ + PHNodeIterator iter(topNode); + PHCompositeNode* dstNode = + dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); + + if (!dstNode) + { + throw std::runtime_error("CombineTowerInfo: DST node not found"); + } + + PHCompositeNode *DetNode = dynamic_cast(iter.findFirst("PHCompositeNode", m_detector)); + + m_towersA = findNode::getClass(topNode, m_inputNodeA); + m_towersB = findNode::getClass(topNode, m_inputNodeB); + + if (!m_towersB) + { + std::cout << "CombineTowerInfo: " <(dstNode, m_outputNode); + if (!m_towersOut) + { + m_towersOut = + dynamic_cast(m_towersA->CloneMe()); + + auto* node = new PHIODataNode( + m_towersOut, m_outputNode, "PHObject"); + + DetNode->addNode(node); + } + + if (m_towersA->size() != m_towersB->size()) + { + throw std::runtime_error("CombineTowerInfo: input containers have different sizes"); + } +} + +//____________________________________________________________________________ +int CombineTowerInfo::process_event(PHCompositeNode* /*topNode*/) +{ + const unsigned int ntowers = m_towersA->size(); + + for (unsigned int ich = 0; ich < ntowers; ++ich) + { + TowerInfo* towerA = m_towersA->get_tower_at_channel(ich); + TowerInfo* towerB = m_towersB->get_tower_at_channel(ich); + TowerInfo* towerO = m_towersOut->get_tower_at_channel(ich); + + towerO->copy_tower(towerA); + + const float e_sum = towerA->get_energy() + towerB->get_energy(); + towerO->set_energy(e_sum); + } + + return Fun4AllReturnCodes::EVENT_OK; +} + diff --git a/offline/packages/CaloEmbedding/CombineTowerInfo.h b/offline/packages/CaloEmbedding/CombineTowerInfo.h new file mode 100644 index 0000000000..ab96a918e6 --- /dev/null +++ b/offline/packages/CaloEmbedding/CombineTowerInfo.h @@ -0,0 +1,39 @@ +#ifndef COMBINETOWERINFO_H +#define COMBINETOWERINFO_H + +#include + +#include + +class PHCompositeNode; +class TowerInfoContainer; + +class CombineTowerInfo : public SubsysReco +{ + public: + explicit CombineTowerInfo(const std::string& name = "CombineTowerInfo"); + ~CombineTowerInfo() override = default; + + int InitRun(PHCompositeNode* topNode) override; + int process_event(PHCompositeNode* topNode) override; + + void set_inputNodeA(const std::string& name) { m_inputNodeA = name; } + void set_inputNodeB(const std::string& name) { m_inputNodeB = name; } + void set_outputNode(const std::string& name) { m_outputNode = name; } + void set_detector(const std::string& name) { m_detector = name; } + + private: + void CreateNodes(PHCompositeNode* topNode); + + std::string m_inputNodeA; + std::string m_inputNodeB; + std::string m_outputNode; + std::string m_detector; + + TowerInfoContainer* m_towersA{nullptr}; + TowerInfoContainer* m_towersB{nullptr}; + TowerInfoContainer* m_towersOut{nullptr}; +}; + +#endif + diff --git a/offline/packages/CaloEmbedding/CopyIODataNodes.cc b/offline/packages/CaloEmbedding/CopyIODataNodes.cc index 0d4d4536ff..08428dccfe 100644 --- a/offline/packages/CaloEmbedding/CopyIODataNodes.cc +++ b/offline/packages/CaloEmbedding/CopyIODataNodes.cc @@ -7,6 +7,9 @@ #include +#include +#include + #include #include @@ -57,6 +60,10 @@ int CopyIODataNodes::InitRun(PHCompositeNode *topNode) { CreateSyncObject(topNode, se->topNode()); } + if (m_CopyTowerInfoFlag) + { + CreateTowerInfo(topNode, se->topNode()); + } return Fun4AllReturnCodes::EVENT_OK; } @@ -89,6 +96,10 @@ int CopyIODataNodes::process_event(PHCompositeNode *topNode) { CopySyncObject(topNode, se->topNode()); } + if (m_CopyTowerInfoFlag) + { + CopyTowerInfo(topNode, se->topNode()); + } return Fun4AllReturnCodes::EVENT_OK; } @@ -293,6 +304,29 @@ void CopyIODataNodes::CopyMinimumBiasInfo(PHCompositeNode *from_topNode, PHCompo return; } +void CopyIODataNodes::CopyTowerInfo(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode) +{ + TowerInfoContainer *from_towerInfo = findNode::getClass(from_topNode, from_towerInfo_name); + TowerInfoContainer *to_towerInfo = findNode::getClass( to_topNode, to_towerInfo_name); + unsigned int ntowers = from_towerInfo->size(); + for (unsigned int ch = 0; ch < ntowers; ++ch) + { + TowerInfo *from_tow = from_towerInfo->get_tower_at_channel(ch); + to_towerInfo->get_tower_at_channel(ch)->copy_tower(from_tow); + } + + if (Verbosity() > 0) + { + std::cout << "From TowerInfoContainer identify()" << std::endl; + from_towerInfo->identify(); + std::cout << "To TowerInfoCOntainer identify()" << std::endl; + to_towerInfo->identify(); + } + + return; +} + + void CopyIODataNodes::CreateMbdOut(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode) { @@ -330,6 +364,36 @@ void CopyIODataNodes::CreateMbdOut(PHCompositeNode *from_topNode, PHCompositeNod } +void CopyIODataNodes::CreateTowerInfo(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode) +{ + std::cout << "copying tower info" << std::endl; + TowerInfoContainer *from_towerInfo = findNode::getClass(from_topNode, from_towerInfo_name); + if (!from_towerInfo) + { + std::cout << "Could not locate TowerInfoContainer on " << from_topNode->getName() << std::endl; + m_CopyTowerInfoFlag = false; + return; + } + TowerInfoContainer *to_towerInfo = findNode::getClass(to_topNode, to_towerInfo_name); + if (!to_towerInfo) + { + PHNodeIterator iter(to_topNode); + PHCompositeNode *dstNode = dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); + if (!dstNode) + { + dstNode = new PHCompositeNode("DST"); + to_topNode->addNode(dstNode); + } + to_towerInfo = dynamic_cast(from_towerInfo->CloneMe()); + PHIODataNode *newNode = new PHIODataNode(to_towerInfo, to_towerInfo_name, "PHObject"); + dstNode->addNode(newNode); + } + return; +} + + + + void CopyIODataNodes::CopyMbdOut(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode) { MbdOut *from_mbdout = findNode::getClass(from_topNode, "MbdOut"); diff --git a/offline/packages/CaloEmbedding/CopyIODataNodes.h b/offline/packages/CaloEmbedding/CopyIODataNodes.h index 08626bdafb..5045792a2d 100644 --- a/offline/packages/CaloEmbedding/CopyIODataNodes.h +++ b/offline/packages/CaloEmbedding/CopyIODataNodes.h @@ -35,6 +35,13 @@ class CopyIODataNodes : public SubsysReco void CopyMbdOut(bool flag = true) { m_CopyMbdOutFlag = flag; } void CopyRunHeader(bool flag = true) { m_CopyRunHeaderFlag = flag; } void CopySyncObject(bool flag = true) { m_CopySyncObjectFlag = flag; } + void set_CopyTowerInfo(const std::string& set_from_towerInfo_name,const std::string& set_to_towerInfo_name) + { + from_towerInfo_name = set_from_towerInfo_name; + to_towerInfo_name = set_to_towerInfo_name; + m_CopyTowerInfoFlag = true; + return; + } private: void CreateCentralityInfo(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode); @@ -56,6 +63,8 @@ class CopyIODataNodes : public SubsysReco void CreateSyncObject(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode); void CopySyncObject(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode); + void CopyTowerInfo(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode); + void CreateTowerInfo(PHCompositeNode *from_topNode, PHCompositeNode *to_topNode); bool m_CopyCentralityInfoFlag = true; bool m_CopyEventHeaderFlag = true; @@ -64,6 +73,10 @@ class CopyIODataNodes : public SubsysReco bool m_CopyMbdOutFlag = true; bool m_CopyRunHeaderFlag = true; bool m_CopySyncObjectFlag = true; + bool m_CopyTowerInfoFlag = false; + + std::string from_towerInfo_name = {}; + std::string to_towerInfo_name = {}; }; #endif // COPYIODATANODES_H diff --git a/offline/packages/CaloEmbedding/Makefile.am b/offline/packages/CaloEmbedding/Makefile.am index 9731b33cb8..f9dd1faa22 100644 --- a/offline/packages/CaloEmbedding/Makefile.am +++ b/offline/packages/CaloEmbedding/Makefile.am @@ -13,6 +13,7 @@ AM_LDFLAGS = \ pkginclude_HEADERS = \ caloTowerEmbed.h \ CopyIODataNodes.h \ + CombineTowerInfo.h \ HepMCCollisionVertex.h lib_LTLIBRARIES = \ @@ -21,6 +22,7 @@ lib_LTLIBRARIES = \ libCaloEmbedding_la_SOURCES = \ caloTowerEmbed.cc \ CopyIODataNodes.cc \ + CombineTowerInfo.cc \ HepMCCollisionVertex.cc libCaloEmbedding_la_LIBADD = \ diff --git a/offline/packages/CaloReco/PhotonClusterBuilder.cc b/offline/packages/CaloReco/PhotonClusterBuilder.cc index 263ae04ef0..edcf5f4418 100644 --- a/offline/packages/CaloReco/PhotonClusterBuilder.cc +++ b/offline/packages/CaloReco/PhotonClusterBuilder.cc @@ -566,6 +566,7 @@ void PhotonClusterBuilder::calculate_shower_shapes(RawCluster* rc, PhotonCluster photon->set_shower_shape_parameter("et3", showershape[2]); photon->set_shower_shape_parameter("et4", showershape[3]); photon->set_shower_shape_parameter("e11", e11); + photon->set_shower_shape_parameter("e22", showershape[8] + showershape[9] + showershape[10] + showershape[11]); photon->set_shower_shape_parameter("e33", e33); photon->set_shower_shape_parameter("e55", e55); photon->set_shower_shape_parameter("e77", e77); diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc index 38b60d55ba..8c049f30d6 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc @@ -148,7 +148,7 @@ int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) { std::cout << "KFParticle: Event skipped as there are no tracks" << std::endl; } - return Fun4AllReturnCodes::ABORTEVENT; + return Fun4AllReturnCodes::EVENT_OK; } if (!m_use_fake_pv) @@ -162,7 +162,7 @@ int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) { std::cout << "KFParticle: Event skipped as there are no vertices" << std::endl; } - return Fun4AllReturnCodes::ABORTEVENT; + return Fun4AllReturnCodes::EVENT_OK; } } else @@ -174,10 +174,9 @@ int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) { std::cout << "KFParticle: Event skipped as there are no vertices" << std::endl; } - return Fun4AllReturnCodes::ABORTEVENT; + return Fun4AllReturnCodes::EVENT_OK; } } - } createDecay(topNode, mother, vertex_kfparticle, daughters, intermediates, nPVs); diff --git a/offline/packages/PHGenFitPkg/PHGenFit/Fitter.cc b/offline/packages/PHGenFitPkg/PHGenFit/Fitter.cc index 8fe89b5019..aeb01dee62 100644 --- a/offline/packages/PHGenFitPkg/PHGenFit/Fitter.cc +++ b/offline/packages/PHGenFitPkg/PHGenFit/Fitter.cc @@ -51,9 +51,9 @@ namespace PHGenFit const std::string& /*track_rep_choice*/, const bool doEventDisplay) : verbosity(1000) + , _tgeo_manager(new TGeoManager("Default", "Geane geometry")) , _doEventDisplay(doEventDisplay) { - _tgeo_manager = new TGeoManager("Default", "Geane geometry"); TGeoManager::Import(tgeo_file_name.data()); assert(field); @@ -73,26 +73,24 @@ namespace PHGenFit } // init fitter - if (fitter_choice.compare("KalmanFitterRefTrack") == 0) + if (fitter_choice == "KalmanFitterRefTrack") { _fitter = new genfit::KalmanFitterRefTrack(); } - else if (fitter_choice.compare("KalmanFitter") == 0) -// NOLINTNEXTLINE(bugprone-branch-clone) - { + else if (fitter_choice == "KalmanFitter") + { // NOLINT(bugprone-branch-clone) _fitter = new genfit::KalmanFitter(); } - else if (fitter_choice.compare("DafSimple") == 0) + else if (fitter_choice == "DafSimple") { _fitter = new genfit::DAF(false); } - else if (fitter_choice.compare("DafRef") == 0) + else if (fitter_choice == "DafRef") { _fitter = new genfit::DAF(true); } else -// NOLINTNEXTLINE(bugprone-branch-clone) - { + { // NOLINT(bugprone-branch-clone) _fitter = new genfit::KalmanFitter(); } @@ -235,7 +233,7 @@ namespace PHGenFit { _fitter = new genfit::KalmanFitterRefTrack(); } - if (fitter_choice == PHGenFit::Fitter::DafSimple) + else if (fitter_choice == PHGenFit::Fitter::DafSimple) { _fitter = new genfit::DAF(false); } @@ -289,19 +287,19 @@ namespace PHGenFit } // init fitter - if (fitter_choice.compare("KalmanFitterRefTrack") == 0) + if (fitter_choice == "KalmanFitterRefTrack") { _fitter = new genfit::KalmanFitterRefTrack(); } - else if (fitter_choice.compare("KalmanFitter") == 0) + else if (fitter_choice == "KalmanFitter") { _fitter = new genfit::KalmanFitter(); } - else if (fitter_choice.compare("DafSimple") == 0) + else if (fitter_choice == "DafSimple") { _fitter = new genfit::DAF(false); } - else if (fitter_choice.compare("DafRef") == 0) + else if (fitter_choice == "DafRef") { _fitter = new genfit::DAF(true); } diff --git a/offline/packages/PHGenFitPkg/PHGenFit/Track.cc b/offline/packages/PHGenFitPkg/PHGenFit/Track.cc index b1643e8cd0..0e06f15984 100644 --- a/offline/packages/PHGenFitPkg/PHGenFit/Track.cc +++ b/offline/packages/PHGenFitPkg/PHGenFit/Track.cc @@ -48,8 +48,8 @@ #define WILD_DOUBLE (-999999) -//#define _DEBUG_ -//#define _PRINT_MATRIX_ +// #define _DEBUG_ +// #define _PRINT_MATRIX_ #ifdef _DEBUG_ #include @@ -60,11 +60,10 @@ ofstream fout_matrix("matrix.txt"); namespace PHGenFit { Track::Track(genfit::AbsTrackRep* rep, const TVector3& seed_pos, const TVector3& seed_mom, const TMatrixDSym& seed_cov, const int v) + : verbosity(v) { // TODO Add input param check - verbosity = v; - genfit::MeasuredStateOnPlane seedMSoP(rep); seedMSoP.setPosMomCov(seed_pos, seed_mom, seed_cov); // const genfit::StateOnPlane seedSoP(seedMSoP); @@ -78,12 +77,12 @@ namespace PHGenFit } Track::Track(const PHGenFit::Track& t) + : verbosity(t.verbosity) + , _track(new genfit::Track(*(t.getGenFitTrack()))) + , _clusterIDs(t.get_cluster_IDs()) + , _clusterkeys(t.get_cluster_keys()) + , _vertex_id(t.get_vertex_id()) { - _track = new genfit::Track(*(t.getGenFitTrack())); - verbosity = t.verbosity; - _clusterIDs = t.get_cluster_IDs(); - _clusterkeys = t.get_cluster_keys(); - _vertex_id = t.get_vertex_id(); } int Track::addMeasurement(PHGenFit::Measurement* measurement) @@ -191,10 +190,8 @@ namespace PHGenFit delete state; return nullptr; } - else - { - return state; - } + + return state; } double Track::extrapolateToLine(genfit::MeasuredStateOnPlane& state, const TVector3& line_point, const TVector3& line_direction, const int tr_point_id) const @@ -240,10 +237,8 @@ namespace PHGenFit delete state; return nullptr; } - else - { - return state; - } + + return state; } double Track::extrapolateToCylinder(genfit::MeasuredStateOnPlane& state, double radius, const TVector3& line_point, const TVector3& line_direction, const int tr_point_id, const int direction) const @@ -361,10 +356,8 @@ namespace PHGenFit delete state; return nullptr; } - else - { - return state; - } + + return state; } int Track::updateOneMeasurementKalman( @@ -385,7 +378,7 @@ namespace PHGenFit << std::endl; #endif - if (measurements.size() == 0) + if (measurements.empty()) { return -1; } @@ -437,11 +430,11 @@ namespace PHGenFit #endif continue; } - //#ifdef _DEBUG_ + // #ifdef _DEBUG_ // std::cout << __LINE__ << "\n ###################################################################"<Print(); // std::cout << __LINE__ << "\n ###################################################################"<getFittedState(true)); @@ -579,7 +572,7 @@ namespace PHGenFit // std::cout << err_phi << "\t" << err_z << "\t"; } #endif - for (auto rawMeasurement : rawMeasurements) + for (auto* rawMeasurement : rawMeasurements) { fi->addMeasurementsOnPlane( rawMeasurement->constructMeasurementsOnPlane(*state)); @@ -598,7 +591,7 @@ namespace PHGenFit << ": size of fi's MeasurementsOnPlane: " << measurements_on_plane.size() << std::endl; #endif - for (auto it : measurements_on_plane) + for (auto* it : measurements_on_plane) { const genfit::MeasurementOnPlane& mOnPlane = *it; // const double weight = mOnPlane.getWeight(); @@ -769,10 +762,7 @@ namespace PHGenFit delete state; return nullptr; } - else - { - return state; - } + return state; } double Track::get_chi2() const diff --git a/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.cc b/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.cc index b519fd472e..224d172a24 100644 --- a/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.cc +++ b/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.cc @@ -20,6 +20,12 @@ TriggerDSTSkimmer::TriggerDSTSkimmer(const std::string &name) //____________________________________________________________________________.. int TriggerDSTSkimmer::process_event(PHCompositeNode *topNode) { + + if ((accepted_events >= max_accept) && use_max_accept) + { + return Fun4AllReturnCodes::ABORTEVENT; + } + if (Verbosity() > 0) { if (ievent % 1000 == 0) @@ -45,7 +51,7 @@ int TriggerDSTSkimmer::process_event(PHCompositeNode *topNode) if (n_trigger_index != 0) { bool trigger_fired = false; - Gl1Packet *_gl1PacketInfo = findNode::getClass(topNode, "GL1Packet"); + Gl1Packet *_gl1PacketInfo = findNode::getClass(topNode, 14001); int gl1_trigger_vector_scaled[64] = {0}; if (_gl1PacketInfo) { @@ -61,6 +67,7 @@ int TriggerDSTSkimmer::process_event(PHCompositeNode *topNode) std::cout << "TriggerDSTSkimmer::process_event - Error - Can't find Trigger Node Gl1Packet therefore no selection can be made" << std::endl; return Fun4AllReturnCodes::ABORTEVENT; } + for (int it = 0; it < n_trigger_index; ++it) { if (gl1_trigger_vector_scaled[m_trigger_index[it]] == 1) @@ -74,5 +81,8 @@ int TriggerDSTSkimmer::process_event(PHCompositeNode *topNode) return Fun4AllReturnCodes::ABORTEVENT; } } + + accepted_events++; + return Fun4AllReturnCodes::EVENT_OK; } diff --git a/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.h b/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.h index dfb2c47a7c..9919df06b7 100644 --- a/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.h +++ b/offline/packages/Skimmers/Trigger/TriggerDSTSkimmer.h @@ -22,10 +22,22 @@ class TriggerDSTSkimmer : public SubsysReco void SetTrigger(std::vector &trigger_vector) {m_trigger_index = trigger_vector;} + void set_accept_max(int max_events) + { + use_max_accept = true; + max_accept = max_events; + return; + } + private: std::vector m_trigger_index{10}; int ievent{0}; + + int accepted_events{0}; + int max_accept{0}; + bool use_max_accept{false}; + }; #endif // JETDSTSKIMMER_H diff --git a/offline/packages/jetbackground/DetermineTowerBackground.cc b/offline/packages/jetbackground/DetermineTowerBackground.cc index 9823d337f6..c359b6914c 100644 --- a/offline/packages/jetbackground/DetermineTowerBackground.cc +++ b/offline/packages/jetbackground/DetermineTowerBackground.cc @@ -14,6 +14,10 @@ #include #include +#include +#include +#include + #include #include @@ -52,9 +56,57 @@ DetermineTowerBackground::DetermineTowerBackground(const std::string &name) int DetermineTowerBackground::InitRun(PHCompositeNode *topNode) { + if (_do_flow == 4) + { + if (Verbosity()) + { + std::cout << "Loading the average calo v2" << std::endl; + } + if (LoadCalibrations()) + { + std::cout << "Load calibrations failed." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + } + return CreateNode(topNode); } +int DetermineTowerBackground::LoadCalibrations() +{ + + CDBTTree *cdbtree_calo_v2 = nullptr; + + std::string calibdir = CDBInterface::instance()->getUrl(m_calibName); + if (m_overwrite_average_calo_v2) + { + calibdir = m_overwrite_average_calo_v2_path; + } + + if (calibdir.empty()) + { + std::cout << "Could not find filename for calo average v2, exiting" << std::endl; + exit(-1); + } + + cdbtree_calo_v2 = new CDBTTree(calibdir); + + + cdbtree_calo_v2->LoadCalibrations(); + + _CENTRALITY_V2.assign(100,0); + + for (int icent = 0; icent < 100; icent++) + { + _CENTRALITY_V2[icent] = cdbtree_calo_v2->GetFloatValue(icent, "jet_calo_v2"); + } + + delete cdbtree_calo_v2; + + return Fun4AllReturnCodes::EVENT_OK; +} + int DetermineTowerBackground::process_event(PHCompositeNode *topNode) { @@ -481,7 +533,92 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) } } - if ( _do_flow >= 1 ) + + // Get psi + if (_do_flow == 2) + { // HIJING truth flow extraction + PHG4TruthInfoContainer *truthinfo = findNode::getClass(topNode, "G4TruthInfo"); + + if (!truthinfo) + { + std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl; + return -1; + } + + PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange(); + + float Hijing_Qx = 0; + float Hijing_Qy = 0; + + for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) + { + PHG4Particle *g4particle = iter->second; + + if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) + { + continue; + } + + TLorentzVector t; + t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e()); + + float truth_pt = t.Pt(); + if (truth_pt < 0.4) + { + continue; + } + float truth_eta = t.Eta(); + if (std::fabs(truth_eta) > 1.1) + { + continue; + } + float truth_phi = t.Phi(); + int truth_pid = g4particle->get_pid(); + + if (Verbosity() > 10) + { + std::cout << "DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt << " / " << truth_eta << " / " << truth_phi << " , embed / PID = " << truthinfo->isEmbeded(g4particle->get_track_id()) << " / " << truth_pid << std::endl; + } + + Hijing_Qx += truth_pt * std::cos(2 * truth_phi); + Hijing_Qy += truth_pt * std::sin(2 * truth_phi); + } + + _Psi2 = std::atan2(Hijing_Qy, Hijing_Qx) / 2.0; + + if (Verbosity() > 0) + { + std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; + } + } + else if (_do_flow == 3 || _do_flow == 4) + { // sEPD event plane extraction + // get event plane map + EventplaneinfoMap *epmap = findNode::getClass(topNode, "EventplaneinfoMap"); + if (!epmap) + { + std::cout << "DetermineTowerBackground::process_event: FATAL, EventplaneinfoMap does not exist, cannot extract sEPD flow with do_flow = " << _do_flow << std::endl; + exit(-1); + } + if (!(epmap->empty())) + { + auto *EPDNS = epmap->get(EventplaneinfoMap::sEPDNS); + _Psi2 = EPDNS->get_shifted_psi(2); + } + else + { + _is_flow_failure = true; + _Psi2 = 0; + } + + if (Verbosity() > 0) + { + std::cout << "DetermineTowerBackground::process_event: flow extracted from sEPD, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; + } + + } + + if ( _do_flow >= 1 && _do_flow < 4) { if (Verbosity() > 0) @@ -754,88 +891,6 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) { // Calo event plane _Psi2 = std::atan2(Q_y, Q_x) / 2.0; } - else if (_do_flow == 2) - { // HIJING truth flow extraction - PHG4TruthInfoContainer *truthinfo = findNode::getClass(topNode, "G4TruthInfo"); - - if (!truthinfo) - { - std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl; - return -1; - } - - PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange(); - - float Hijing_Qx = 0; - float Hijing_Qy = 0; - - for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) - { - PHG4Particle *g4particle = iter->second; - - if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) - { - continue; - } - - TLorentzVector t; - t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e()); - - float truth_pt = t.Pt(); - if (truth_pt < 0.4) - { - continue; - } - float truth_eta = t.Eta(); - if (std::fabs(truth_eta) > 1.1) - { - continue; - } - float truth_phi = t.Phi(); - int truth_pid = g4particle->get_pid(); - - if (Verbosity() > 10) - { - std::cout << "DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt << " / " << truth_eta << " / " << truth_phi << " , embed / PID = " << truthinfo->isEmbeded(g4particle->get_track_id()) << " / " << truth_pid << std::endl; - } - - Hijing_Qx += truth_pt * std::cos(2 * truth_phi); - Hijing_Qy += truth_pt * std::sin(2 * truth_phi); - } - - _Psi2 = std::atan2(Hijing_Qy, Hijing_Qx) / 2.0; - - if (Verbosity() > 0) - { - std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; - } - } - else if (_do_flow == 3) - { // sEPD event plane extraction - // get event plane map - EventplaneinfoMap *epmap = findNode::getClass(topNode, "EventplaneinfoMap"); - if (!epmap) - { - std::cout << "DetermineTowerBackground::process_event: FATAL, EventplaneinfoMap does not exist, cannot extract sEPD flow with do_flow = " << _do_flow << std::endl; - exit(-1); - } - if (!(epmap->empty())) - { - auto *EPDNS = epmap->get(EventplaneinfoMap::sEPDNS); - _Psi2 = EPDNS->get_shifted_psi(2); - } - else - { - _is_flow_failure = true; - _Psi2 = 0; - } - - if (Verbosity() > 0) - { - std::cout << "DetermineTowerBackground::process_event: flow extracted from sEPD, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; - } - - } if (std::isnan(_Psi2) || std::isinf(_Psi2)) { @@ -890,7 +945,30 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) std::cout << "DetermineTowerBackground::process_event: flow extraction successful, Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) , v2 = " << _v2 << std::endl; } } // if do flow + else if (_do_flow == 4) + { + CentralityInfo *centinfo = findNode::getClass(topNode, "CentralityInfo"); + + if (!centinfo) + { + std::cout << "DetermineTowerBackground::process_event: FATAL, CentralityInfo does not exist, cannot extract centrality with do_flow = " << _do_flow << std::endl; + exit(-1); + } + + int centrality_bin = centinfo->get_centrality_bin(CentralityInfo::PROP::mbd_NS); + + if (centrality_bin > 0 && centrality_bin < 95) + { + _v2 = _CENTRALITY_V2[centrality_bin]; + } + else + { + _v2 = 0; + _is_flow_failure = true; + _Psi2 = 0; + } + } // now calculate energy densities... _nTowers = 0; // store how many towers were used to determine bkg diff --git a/offline/packages/jetbackground/DetermineTowerBackground.h b/offline/packages/jetbackground/DetermineTowerBackground.h index a8a6d0209c..ed9e34bf9a 100644 --- a/offline/packages/jetbackground/DetermineTowerBackground.h +++ b/offline/packages/jetbackground/DetermineTowerBackground.h @@ -13,6 +13,7 @@ #include #include #include +#include // forward declarations class PHCompositeNode; @@ -37,7 +38,11 @@ class DetermineTowerBackground : public SubsysReco void SetBackgroundOutputName(const std::string &name) { _backgroundName = name; } void SetSeedType(int seed_type) { _seed_type = seed_type; } void SetFlow(int do_flow) { _do_flow = do_flow; }; - + void SetOverwriteCaloV2(std::string &url) + { + m_overwrite_average_calo_v2 = true; + m_overwrite_average_calo_v2_path = url; + } void SetSeedJetD(float D) { _seed_jet_D = D; }; void SetSeedJetPt(float pt) { _seed_jet_pt = pt; }; void SetSeedMaxConst(float max_const) { _seed_max_const = max_const; }; @@ -55,6 +60,13 @@ class DetermineTowerBackground : public SubsysReco int CreateNode(PHCompositeNode *topNode); void FillNode(PHCompositeNode *topNode); + int LoadCalibrations(); + + std::vector _CENTRALITY_V2; + std::string m_calibName = "JET_AVERAGE_CALO_V2_SEPD_PSI2"; + bool m_overwrite_average_calo_v2{false}; + std::string m_overwrite_average_calo_v2_path; + int _do_flow{0}; float _v2{0}; float _Psi2{0}; diff --git a/offline/packages/jetbackground/Makefile.am b/offline/packages/jetbackground/Makefile.am index 6fe9740a4e..ad45676fbe 100644 --- a/offline/packages/jetbackground/Makefile.am +++ b/offline/packages/jetbackground/Makefile.am @@ -24,6 +24,8 @@ libjetbackground_la_LDFLAGS = \ libjetbackground_la_LIBADD = \ libjetbackground_io.la \ -lcalo_io \ + -lcentrality_io \ + -lcdbobjects \ -lConstituentSubtractor \ -leventplaneinfo_io \ -lglobalvertex \ @@ -33,6 +35,7 @@ libjetbackground_la_LIBADD = \ -lphg4hit \ -lphparameter \ -lqautils \ + -lffamodules \ -lSubsysReco pkginclude_HEADERS = \ diff --git a/offline/packages/mbd/MbdCalib.cc b/offline/packages/mbd/MbdCalib.cc index 22b1e9930d..0b02908adb 100644 --- a/offline/packages/mbd/MbdCalib.cc +++ b/offline/packages/mbd/MbdCalib.cc @@ -1331,7 +1331,7 @@ int MbdCalib::Download_TimeRMS(const std::string& dbase_location) if ( _trms_y[0].empty() ) { - std::cout << PHWHERE << ", ERROR, unknown file type, " << dbase_location << std::endl; + std::cout << PHWHERE << ", WARNING, trms calib missing " << dbase_location << std::endl; _status = -1; return _status; // file not found } diff --git a/offline/packages/mbd/MbdEvent.cc b/offline/packages/mbd/MbdEvent.cc index d8c72b4395..cdab537f45 100644 --- a/offline/packages/mbd/MbdEvent.cc +++ b/offline/packages/mbd/MbdEvent.cc @@ -451,6 +451,7 @@ int MbdEvent::SetRawData(std::array< CaloPacket *,2> &dstp, MbdRawContainer *bbc if (dstp[ipkt]) { _nsamples = dstp[ipkt]->iValue(0, "SAMPLES"); + { static bool printcount{true}; if ( printcount && Verbosity() > 0) @@ -460,6 +461,13 @@ int MbdEvent::SetRawData(std::array< CaloPacket *,2> &dstp, MbdRawContainer *bbc } } + // skip empty packets, corrupt event + if ( _nsamples == 0 ) + { + std::cout << PHWHERE << " ERROR, evt " << m_evt << " no samples in Packet " << pktid << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + m_xmitclocks[ipkt] = static_cast(dstp[ipkt]->iValue(0, "CLOCK")); m_femclocks[ipkt][0] = static_cast(dstp[ipkt]->iValue(0, "FEMCLOCK")); @@ -484,9 +492,17 @@ int MbdEvent::SetRawData(std::array< CaloPacket *,2> &dstp, MbdRawContainer *bbc } _mbdsig[feech].SetNSamples( _nsamples ); - _mbdsig[feech].SetXY(m_samp[feech], m_adc[feech]); - + + if ( _nsamples > 0 && _nsamples <= 30 ) + { + _mbdsig[feech].SetXY(m_samp[feech], m_adc[feech]); + } /* + else + { + std::cout << PHWHERE << " empty feech " << feech << std::endl; + } + std::cout << "feech " << feech << std::endl; _mbdsig[feech].Print(); */ @@ -565,6 +581,7 @@ int MbdEvent::SetRawData(Event *event, MbdRawContainer *bbcraws, MbdPmtContainer if (p[ipkt]) { _nsamples = p[ipkt]->iValue(0, "SAMPLES"); + { static int counter = 0; if ( counter<1 ) @@ -574,6 +591,15 @@ int MbdEvent::SetRawData(Event *event, MbdRawContainer *bbcraws, MbdPmtContainer counter++; } + // If packets are missing, stop processing event + if ( _nsamples == 0 ) + { + std::cout << PHWHERE << " ERROR, skipping evt " << m_evt << " nsamples = 0 " << pktid << std::endl; + delete p[ipkt]; + p[ipkt] = nullptr; + return Fun4AllReturnCodes::ABORTEVENT; + } + m_xmitclocks[ipkt] = static_cast(p[ipkt]->iValue(0, "CLOCK")); m_femclocks[ipkt][0] = static_cast(p[ipkt]->iValue(0, "FEMCLOCK")); @@ -637,7 +663,7 @@ int MbdEvent::ProcessPackets(MbdRawContainer *bbcraws) // Do a quick sanity check that all fem counters agree if (m_xmitclocks[0] != m_xmitclocks[1]) { - std::cout << __FILE__ << ":" << __LINE__ << " ERROR, xmitclocks don't agree" << std::endl; + std::cout << __FILE__ << ":" << __LINE__ << " ERROR, xmitclocks don't agree, evt " << m_evt << std::endl; } /* // format changed in run2024, need to update check @@ -673,6 +699,11 @@ int MbdEvent::ProcessPackets(MbdRawContainer *bbcraws) int pmtch = _mbdgeom->get_pmt(ifeech); int type = _mbdgeom->get_type(ifeech); // 0 = T-channel, 1 = Q-channel + if ( _mbdsig[ifeech].GetNSamples()==0 ) + { + continue; + } + // time channel if (type == 0) { @@ -739,6 +770,11 @@ int MbdEvent::ProcessRawContainer(MbdRawContainer *bbcraws, MbdPmtContainer *bbc int pmtch = _mbdgeom->get_pmt(ifeech); int type = _mbdgeom->get_type(ifeech); // 0 = T-channel, 1 = Q-channel + if ( _mbdsig[ifeech].GetNSamples()==0 ) + { + continue; + } + // time channel if (type == 0) { @@ -854,8 +890,14 @@ int MbdEvent::ProcessRawContainer(MbdRawContainer *bbcraws, MbdPmtContainer *bbc */ TGraphErrors *gsubpulse = _mbdsig[ifeech].GetGraph(); - Double_t *y = gsubpulse->GetY(); - h2_trange->Fill( y[samp_max], pmtch ); // fill ped-subtracted tdc + if ( gsubpulse ) + { + Double_t *y = gsubpulse->GetY(); + if ( y ) + { + h2_trange->Fill( y[samp_max], pmtch ); // fill ped-subtracted tdc + } + } } } diff --git a/offline/packages/mbd/MbdReco.cc b/offline/packages/mbd/MbdReco.cc index c19ed37ba9..004e191e78 100644 --- a/offline/packages/mbd/MbdReco.cc +++ b/offline/packages/mbd/MbdReco.cc @@ -103,7 +103,8 @@ int MbdReco::process_event(PHCompositeNode *topNode) int status = Fun4AllReturnCodes::EVENT_OK; if ( m_evtheader!=nullptr ) { - m_mbdevent->set_EventNumber( m_evtheader->get_EvtSequence() ); + _evtnum = m_evtheader->get_EvtSequence(); + m_mbdevent->set_EventNumber( _evtnum ); } if ( m_event!=nullptr ) @@ -125,7 +126,7 @@ int MbdReco::process_event(PHCompositeNode *topNode) static int counter = 0; if ( counter<3 ) { - std::cout << PHWHERE << " Warning, MBD discarding event " << std::endl; + std::cout << PHWHERE << " Warning, MBD discarding event " << _evtnum << std::endl; counter++; } return Fun4AllReturnCodes::DISCARDEVENT; @@ -135,7 +136,7 @@ int MbdReco::process_event(PHCompositeNode *topNode) static int counter = 0; if ( counter<3 ) { - std::cout << PHWHERE << " Warning, MBD aborting event " << std::endl; + std::cout << PHWHERE << " Warning, MBD aborting event " << _evtnum << std::endl; counter++; } return Fun4AllReturnCodes::ABORTEVENT; diff --git a/offline/packages/mbd/MbdReco.h b/offline/packages/mbd/MbdReco.h index 0dfc7d7cbc..6e00de68f3 100644 --- a/offline/packages/mbd/MbdReco.h +++ b/offline/packages/mbd/MbdReco.h @@ -52,6 +52,8 @@ class MbdReco : public SubsysReco float m_tres = 0.05; std::unique_ptr m_gaussian = nullptr; + int _evtnum{-1}; + std::unique_ptr m_mbdevent{nullptr}; Event *m_event{nullptr}; std::arraym_mbdpacket{nullptr}; diff --git a/offline/packages/mbd/MbdSig.h b/offline/packages/mbd/MbdSig.h index 0b8b420ffa..e933f33392 100644 --- a/offline/packages/mbd/MbdSig.h +++ b/offline/packages/mbd/MbdSig.h @@ -32,6 +32,8 @@ class MbdSig void SetY(const Float_t *y, const int invert = 1); void SetXY(const Float_t *x, const Float_t *y, const int invert = 1); + int GetNSamples() { return _nsamples; } + void SetCalib(MbdCalib *mcal); TH1 *GetHist() { return hpulse; } diff --git a/offline/packages/micromegas/MicromegasClusterizer.cc b/offline/packages/micromegas/MicromegasClusterizer.cc index c57d6ed3f5..5da9264c11 100644 --- a/offline/packages/micromegas/MicromegasClusterizer.cc +++ b/offline/packages/micromegas/MicromegasClusterizer.cc @@ -157,8 +157,13 @@ int MicromegasClusterizer::process_event(PHCompositeNode *topNode) // geometry PHG4CylinderGeomContainer* geonode = nullptr; for( std::string geonodename: {"CYLINDERGEOM_MICROMEGAS_FULL", "CYLINDERGEOM_MICROMEGAS" } ) - { if(( geonode = findNode::getClass(topNode, geonodename.c_str()) )) { break; -}} + { + // try load node and test + geonode = findNode::getClass(topNode, geonodename); + if( geonode ) { break;} + } + + //ma assert(geonode); // hitset container @@ -182,8 +187,8 @@ int MicromegasClusterizer::process_event(PHCompositeNode *topNode) for( auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it ) { // get hitset, key and layer - TrkrHitSet* hitset = hitset_it->second; - const TrkrDefs::hitsetkey hitsetkey = hitset_it->first; + const auto& [hitsetkey, hitset] = *hitset_it; + const auto layer = TrkrDefs::getLayer(hitsetkey); const auto tileid = MicromegasDefs::getTileId(hitsetkey); @@ -215,17 +220,32 @@ int MicromegasClusterizer::process_event(PHCompositeNode *topNode) using range_list_t = std::vector; range_list_t ranges; - // loop over hits - const auto hit_range = hitset->getHits(); + // Make a local copy of hitsets, sorted along strips + /* when there are multiple hits on the same strip, only the first one (in time) is kept */ + class StripSortFtor + { + public: + bool operator() ( const TrkrDefs::hitkey& first, const TrkrDefs::hitkey& second ) const + { return MicromegasDefs::getStrip(first) < MicromegasDefs::getStrip(second); } + }; + + using LocalMap = std::map; + LocalMap local_hitmap; + + { + // loop over hits + const auto hit_range = hitset->getHits(); + std::copy( hit_range.first, hit_range.second, std::inserter(local_hitmap, local_hitmap.end()) ); + } // keep track of first iterator of runing cluster - auto begin = hit_range.first; + auto begin = local_hitmap.begin(); // keep track of previous strip uint16_t previous_strip = 0; bool first = true; - for( auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it ) + for( auto hit_it = local_hitmap.begin(); hit_it != local_hitmap.end(); ++hit_it ) { // get hit key @@ -233,18 +253,11 @@ int MicromegasClusterizer::process_event(PHCompositeNode *topNode) // get strip number const auto strip = MicromegasDefs::getStrip( hitkey ); - - if( first ) + if( !first && (strip - previous_strip > 1 ) ) { - previous_strip = strip; - first = false; - continue; - - } else if( strip - previous_strip > 1 ) { - // store current cluster range - ranges.push_back( std::make_pair( begin, hit_it ) ); + ranges.emplace_back( begin, hit_it ); // reinitialize begin of next cluster range begin = hit_it; @@ -252,13 +265,13 @@ int MicromegasClusterizer::process_event(PHCompositeNode *topNode) } // update previous strip + first = false; previous_strip = strip; } // store last cluster - if( begin != hit_range.second ) { ranges.push_back( std::make_pair( begin, hit_range.second ) ); -} + if( begin != local_hitmap.end() ) { ranges.emplace_back( begin, local_hitmap.end() ); } // initialize cluster count int cluster_count = 0; diff --git a/offline/packages/micromegas/MicromegasCombinedDataDecoder.cc b/offline/packages/micromegas/MicromegasCombinedDataDecoder.cc index a5269617f7..1c34357280 100644 --- a/offline/packages/micromegas/MicromegasCombinedDataDecoder.cc +++ b/offline/packages/micromegas/MicromegasCombinedDataDecoder.cc @@ -203,13 +203,14 @@ int MicromegasCombinedDataDecoder::process_event(PHCompositeNode* topNode) // loop over sample_range find maximum const auto sample_range = std::make_pair(rawhit->get_sample_begin(), rawhit->get_sample_end()); - std::vector adc_list; + using sample_pair_t = std::pair; + std::vector adc_list; for (auto is = std::max(m_sample_min, sample_range.first); is < std::min(m_sample_max, sample_range.second); ++is) { const uint16_t adc = rawhit->get_adc(is); if (adc != MicromegasDefs::m_adc_invalid) { - adc_list.push_back(adc); + adc_list.emplace_back(is, adc); } } @@ -220,16 +221,18 @@ int MicromegasCombinedDataDecoder::process_event(PHCompositeNode* topNode) // get max adc value in range /* TODO: use more advanced signal processing */ - auto max_adc = *std::max_element(adc_list.begin(), adc_list.end()); + auto max_adc = *std::max_element(adc_list.begin(), adc_list.end(), + [](const sample_pair_t& first, const sample_pair_t& second) + { return first.second < second.second; } ); // compare to hard min_adc value - if (max_adc < m_min_adc) + if (max_adc.second < m_min_adc) { continue; } // compare to threshold - if (max_adc < pedestal + m_n_sigma * rms) + if (max_adc.second < pedestal + m_n_sigma * rms) { continue; } @@ -243,7 +246,8 @@ int MicromegasCombinedDataDecoder::process_event(PHCompositeNode* topNode) << " tile: " << tile << " channel: " << channel << " strip: " << strip - << " adc: " << max_adc + << " sample: " << max_adc.first + << " adc: " << max_adc.second << std::endl; } @@ -251,19 +255,19 @@ int MicromegasCombinedDataDecoder::process_event(PHCompositeNode* topNode) const auto hitset_it = trkrhitsetcontainer->findOrAddHitSet(hitsetkey); // generate hit key - const TrkrDefs::hitkey hitkey = MicromegasDefs::genHitKey(strip); + const TrkrDefs::hitkey hitkey = MicromegasDefs::genHitKey(strip, max_adc.first); // find existing hit, or create - auto hit = hitset_it->second->getHit(hitkey); + auto* hit = hitset_it->second->getHit(hitkey); if (hit) { - // std::cout << "MicromegasCombinedDataDecoder::process_event - duplicated hit, hitsetkey: " << hitsetkey << " strip: " << strip << std::endl; + std::cout << "MicromegasCombinedDataDecoder::process_event - duplicated hit, hitsetkey: " << hitsetkey << " strip: " << strip << std::endl; continue; } // create hit, assign adc and insert in hitset hit = new TrkrHitv2; - hit->setAdc(max_adc); + hit->setAdc(max_adc.second); hitset_it->second->addHitSpecificKey(hitkey, hit); // increment counter diff --git a/offline/packages/micromegas/MicromegasCombinedDataDecoder.h b/offline/packages/micromegas/MicromegasCombinedDataDecoder.h index 32eebc56e0..33033f1992 100644 --- a/offline/packages/micromegas/MicromegasCombinedDataDecoder.h +++ b/offline/packages/micromegas/MicromegasCombinedDataDecoder.h @@ -49,10 +49,10 @@ class MicromegasCombinedDataDecoder : public SubsysReco /** This removes faulty channels for which calibration has failed */ void set_min_adc(double value) { m_min_adc = value; } - /// set min sample for noise estimation + /// set min sample for signal hits void set_sample_min(uint16_t value) { m_sample_min = value; } - /// set min sample for noise estimation + /// set max sample for signal hits void set_sample_max(uint16_t value) { m_sample_max = value; } private: @@ -85,7 +85,7 @@ class MicromegasCombinedDataDecoder : public SubsysReco uint16_t m_sample_min = 0; /// max sample for signal - uint16_t m_sample_max = 100; + uint16_t m_sample_max = 1024; /// keep track of number of hits per hitsetid using hitcountmap_t = std::map; diff --git a/offline/packages/micromegas/MicromegasDefs.cc b/offline/packages/micromegas/MicromegasDefs.cc index 0766d78b8e..8a35461492 100644 --- a/offline/packages/micromegas/MicromegasDefs.cc +++ b/offline/packages/micromegas/MicromegasDefs.cc @@ -25,11 +25,12 @@ namespace * 8 - 16 segmentation type * 0 - 8 tile id */ - static constexpr unsigned int kBitShiftSegmentation = 8; - static constexpr unsigned int kBitShiftTileId = 0; + constexpr unsigned int kBitShiftSegmentation = 8; + constexpr unsigned int kBitShiftTileId = 0; //! bit shift for hit key - static constexpr unsigned int kBitShiftStrip = 0; + constexpr unsigned int kBitShiftStrip = 0; + constexpr unsigned int kBitShiftSample = 8; } @@ -41,10 +42,10 @@ namespace MicromegasDefs { TrkrDefs::hitsetkey key = TrkrDefs::genHitSetKey(TrkrDefs::TrkrId::micromegasId, layer); - TrkrDefs::hitsetkey tmp = to_underlying_type(type); + TrkrDefs::hitsetkey tmp = to_underlying_type(type)&0x1U; key |= (tmp << kBitShiftSegmentation); - tmp = tile; + tmp = tile&0xFFU; key |= (tmp << kBitShiftTileId); return key; @@ -54,28 +55,36 @@ namespace MicromegasDefs SegmentationType getSegmentationType(TrkrDefs::hitsetkey key) { TrkrDefs::hitsetkey tmp = (key >> kBitShiftSegmentation); - return static_cast(tmp); + return static_cast(tmp&0x1U); } //________________________________________________________________ uint8_t getTileId(TrkrDefs::hitsetkey key) { TrkrDefs::hitsetkey tmp = (key >> kBitShiftTileId); - return tmp; + return tmp&0xFFU; } //________________________________________________________________ - TrkrDefs::hitkey genHitKey(uint16_t strip) + TrkrDefs::hitkey genHitKey(uint16_t strip, uint16_t sample) { - TrkrDefs::hitkey key = strip << kBitShiftStrip; - return key; + const TrkrDefs::hitkey key = (strip&0xFFU) << kBitShiftStrip; + const TrkrDefs::hitkey tmp = (sample&0xFFFFU) << kBitShiftSample; + return key|tmp; } //________________________________________________________________ - uint16_t getStrip( TrkrDefs::hitkey key ) + uint8_t getStrip( TrkrDefs::hitkey key ) { TrkrDefs::hitkey tmp = (key >> kBitShiftStrip); - return tmp; + return tmp & 0xFFU; + } + + //________________________________________________________________ + uint16_t getSample( TrkrDefs::hitkey key ) + { + TrkrDefs::hitkey tmp = (key >> kBitShiftSample); + return tmp & 0xFFFFU; } //________________________________________________________________ diff --git a/offline/packages/micromegas/MicromegasDefs.h b/offline/packages/micromegas/MicromegasDefs.h index c95fdffd72..b206b8f895 100644 --- a/offline/packages/micromegas/MicromegasDefs.h +++ b/offline/packages/micromegas/MicromegasDefs.h @@ -60,11 +60,15 @@ namespace MicromegasDefs /*! * @brief Generate a hitkey from strip index inside tile * @param[in] strip strip index + * @param[in] sample sample index */ - TrkrDefs::hitkey genHitKey(uint16_t strip ); + TrkrDefs::hitkey genHitKey(uint16_t strip, uint16_t sample = 0 ); //! get strip from hit key - uint16_t getStrip(TrkrDefs::hitkey); + uint8_t getStrip(TrkrDefs::hitkey); + + //! get sample from hit key + uint16_t getSample(TrkrDefs::hitkey); /*! * @brief Get the segmentation type from cluster key diff --git a/offline/packages/mvtx/CylinderGeom_Mvtx.cc b/offline/packages/mvtx/CylinderGeom_Mvtx.cc index e220d1474d..ab22c80bdc 100644 --- a/offline/packages/mvtx/CylinderGeom_Mvtx.cc +++ b/offline/packages/mvtx/CylinderGeom_Mvtx.cc @@ -8,7 +8,6 @@ #include #include // for operator<<, basic_ostream::operator<<, basic_... -using namespace std; using Segmentation = SegmentationAlpide; CylinderGeom_Mvtx::CylinderGeom_Mvtx( @@ -89,7 +88,7 @@ void CylinderGeom_Mvtx::get_sensor_indices_from_world_coords(std::vector double chip_delta_z = (inner_loc_chip_in_module[8][2] - inner_loc_chip_in_module[0][2]) / 8.0; // int chip_tmp = (int) (world[2]/chip_delta_z) + 4; // 0-9 int chip_tmp = round(world[2] / chip_delta_z) + 4; // 0-9 - // std::cout << " z " << world[2] << " chip_delta_z " << chip_delta_z << " chip_tmp " << chip_tmp << endl; + // std::cout << " z " << world[2] << " chip_delta_z " << chip_delta_z << " chip_tmp " << chip_tmp << std::endl; stave_index = stave_tmp; chip_index = chip_tmp; @@ -102,15 +101,15 @@ bool CylinderGeom_Mvtx::get_pixel_from_local_coords(TVector3 sensor_local, int& double EPS = 5e-6; if (fabs(fabs(sensor_local.X()) - SegmentationAlpide::ActiveMatrixSizeRows / 2.F) < EPS) { - // cout << " Adjusting X, before X= " << sensor_local.X() << endl; + // std::cout << " Adjusting X, before X= " << sensor_local.X() << std::endl; sensor_local.SetX(((sensor_local.X() < 0) ? -1 : 1) * (SegmentationAlpide::ActiveMatrixSizeRows / 2.F - EPS)); - // cout << " Adjusting X, after X= " << sensor_local.X() << endl; + // std::cout << " Adjusting X, after X= " << sensor_local.X() << std::endl; } if (fabs(fabs(sensor_local.Z()) - SegmentationAlpide::ActiveMatrixSizeCols / 2.F) < EPS) { - // cout << " Adjusting Z, before Z= " << sensor_local.Z() << endl; + // std::cout << " Adjusting Z, before Z= " << sensor_local.Z() << std::endl; sensor_local.SetZ(((sensor_local.Z() < 0) ? -1 : 1) * (SegmentationAlpide::ActiveMatrixSizeCols / 2.F - EPS)); - // cout << " Adjusting Z, after Z= " << sensor_local.Z() << endl; + // std::cout << " Adjusting Z, after Z= " << sensor_local.Z() << std::endl; } // YCM (2020-01-02): go from sensor to chip local coords TVector3 in_chip = sensor_local; @@ -122,21 +121,22 @@ bool CylinderGeom_Mvtx::get_pixel_from_local_coords(TVector3 sensor_local, int& int CylinderGeom_Mvtx::get_pixel_from_local_coords(const TVector3& sensor_local) { - int Ngridx, Ngridz; + int Ngridx; + int Ngridz; bool px_in = get_pixel_from_local_coords(sensor_local, Ngridx, Ngridz); if (!px_in) { - cout << PHWHERE + std::cout << PHWHERE << " Pixel is out sensor. (" << sensor_local.X() << ", " << sensor_local.Y() << ", " << sensor_local.Z() << ")." - << endl; + << std::endl; } if (Ngridx < 0 || Ngridx >= get_NX() || Ngridz < 0 || Ngridz >= get_NZ()) { - cout << PHWHERE << "Wrong pixel value X= " << Ngridx << " and Z= " << Ngridz << endl; + std::cout << PHWHERE << "Wrong pixel value X= " << Ngridx << " and Z= " << Ngridz << std::endl; } // numbering starts at zero @@ -157,8 +157,8 @@ TVector3 CylinderGeom_Mvtx::get_local_coords_from_pixel(int iRow, int iCol) bool check = SegmentationAlpide::detectorToLocal((float) iRow, (float) iCol, local); if (!check) { - cout << PHWHERE << "Pixel coord ( " << iRow << ", " << iCol << " )" - << "out of range" << endl; + std::cout << PHWHERE << "Pixel coord ( " << iRow << ", " << iCol << " )" + << "out of range" << std::endl; } // Transform location in chip to location in sensors TVector3 trChipToSens(loc_sensor_in_chip[0], @@ -177,7 +177,7 @@ void CylinderGeom_Mvtx::identify(std::ostream& os) const << ", pixel_x: " << pixel_x << ", pixel_z: " << pixel_z << ", pixel_thickness: " << pixel_thickness - << endl; + << std::endl; return; } @@ -192,17 +192,17 @@ int CylinderGeom_Mvtx::get_NX() const return SegmentationAlpide::NRows; } -int CylinderGeom_Mvtx::get_pixel_X_from_pixel_number(int NXZ) +int CylinderGeom_Mvtx::get_pixel_X_from_pixel_number(int NXZ) const { return NXZ % get_NX(); } -int CylinderGeom_Mvtx::get_pixel_Z_from_pixel_number(int NXZ) +int CylinderGeom_Mvtx::get_pixel_Z_from_pixel_number(int NXZ) const { return NXZ / get_NX(); } -int CylinderGeom_Mvtx::get_pixel_number_from_xbin_zbin(int xbin, int zbin) // obsolete +int CylinderGeom_Mvtx::get_pixel_number_from_xbin_zbin(int xbin, int zbin) const // obsolete { return xbin + zbin * get_NX(); } diff --git a/offline/packages/mvtx/CylinderGeom_Mvtx.h b/offline/packages/mvtx/CylinderGeom_Mvtx.h index e28da2c594..047673a549 100644 --- a/offline/packages/mvtx/CylinderGeom_Mvtx.h +++ b/offline/packages/mvtx/CylinderGeom_Mvtx.h @@ -13,7 +13,7 @@ class CylinderGeom_Mvtx : public PHG4CylinderGeom public: CylinderGeom_Mvtx( int layer, - int in_Nstaves, + int in_N_staves, double in_layer_nominal_radius, double in_phistep, double in_phitilt, @@ -31,7 +31,7 @@ class CylinderGeom_Mvtx : public PHG4CylinderGeom { } - ~CylinderGeom_Mvtx() override {} + ~CylinderGeom_Mvtx() override = default; // from PHObject void identify(std::ostream& os = std::cout) const override; @@ -53,11 +53,11 @@ class CylinderGeom_Mvtx : public PHG4CylinderGeom TVector3 get_local_coords_from_pixel(int NXZ); TVector3 get_local_coords_from_pixel(int iRow, int iCol); - int get_pixel_X_from_pixel_number(int NXZ); + int get_pixel_X_from_pixel_number(int NXZ) const; - int get_pixel_Z_from_pixel_number(int NXZ); + int get_pixel_Z_from_pixel_number(int NXZ) const; - int get_pixel_number_from_xbin_zbin(int xbin, int zbin); // obsolete + int get_pixel_number_from_xbin_zbin(int xbin, int zbin) const; // obsolete double get_stave_phi_tilt() const { return stave_phi_tilt; } double get_stave_phi_0() const { return stave_phi_0; } diff --git a/offline/packages/mvtx/MvtxClusterPruner.cc b/offline/packages/mvtx/MvtxClusterPruner.cc index 24d41db63d..456b302dd9 100644 --- a/offline/packages/mvtx/MvtxClusterPruner.cc +++ b/offline/packages/mvtx/MvtxClusterPruner.cc @@ -12,108 +12,116 @@ #include #include -#include #include +#include -#include #include +#include namespace { //! range adaptor to be able to use range-based for loop - template class range_adaptor + template + class range_adaptor { - public: - range_adaptor( const T& range ):m_range(range){} - const typename T::first_type& begin() {return m_range.first;} - const typename T::second_type& end() {return m_range.second;} - private: + public: + explicit range_adaptor(const T& range) + : m_range(range) + { + } + const typename T::first_type& begin() { return m_range.first; } + const typename T::second_type& end() { return m_range.second; } + + private: T m_range; }; // print cluster information - void print_cluster_information( TrkrDefs::cluskey ckey, TrkrCluster* cluster ) + void print_cluster_information(TrkrDefs::cluskey ckey, TrkrCluster* cluster) { - if( cluster ) + if (cluster) { std::cout << " MVTX cluster: " << ckey - << " position: (" << cluster->getLocalX() << ", " << cluster->getLocalY() << ")" - << " size: " << (int)cluster->getSize() - << " layer: " << (int)TrkrDefs::getLayer(ckey) - << " stave: " << (int) MvtxDefs::getStaveId(ckey) - << " chip: " << (int)MvtxDefs::getChipId(ckey) - << " strobe: " << (int)MvtxDefs::getStrobeId(ckey) - << " index: " << (int)TrkrDefs::getClusIndex(ckey) - << std::endl; - } else { + << " position: (" << cluster->getLocalX() << ", " << cluster->getLocalY() << ")" + << " size: " << (int) cluster->getSize() + << " layer: " << (int) TrkrDefs::getLayer(ckey) + << " stave: " << (int) MvtxDefs::getStaveId(ckey) + << " chip: " << (int) MvtxDefs::getChipId(ckey) + << " strobe: " << MvtxDefs::getStrobeId(ckey) + << " index: " << (int) TrkrDefs::getClusIndex(ckey) + << std::endl; + } + else + { std::cout << " MVTX cluster: " << ckey - << " layer: " << (int)TrkrDefs::getLayer(ckey) - << " stave: " << (int) MvtxDefs::getStaveId(ckey) - << " chip: " << (int)MvtxDefs::getChipId(ckey) - << " strobe: " << (int)MvtxDefs::getStrobeId(ckey) - << " index: " << (int)TrkrDefs::getClusIndex(ckey) - << std::endl; + << " layer: " << (int) TrkrDefs::getLayer(ckey) + << " stave: " << (int) MvtxDefs::getStaveId(ckey) + << " chip: " << (int) MvtxDefs::getChipId(ckey) + << " strobe: " << MvtxDefs::getStrobeId(ckey) + << " index: " << (int) TrkrDefs::getClusIndex(ckey) + << std::endl; } } using hitkeyset_t = std::set; - using clustermap_t = std::map; + using clustermap_t = std::map; -} +} // namespace //_____________________________________________________________________________ -MvtxClusterPruner::MvtxClusterPruner(const std::string &name) +MvtxClusterPruner::MvtxClusterPruner(const std::string& name) : SubsysReco(name) { } //_____________________________________________________________________________ -int MvtxClusterPruner::InitRun(PHCompositeNode * /*topNode*/) +int MvtxClusterPruner::InitRun(PHCompositeNode* /*topNode*/) { std::cout << "MvtxClusterPruner::InitRun - m_use_strict_matching: " << m_use_strict_matching << std::endl; return Fun4AllReturnCodes::EVENT_OK; } //_____________________________________________________________________________ -int MvtxClusterPruner::process_event(PHCompositeNode *topNode) +int MvtxClusterPruner::process_event(PHCompositeNode* topNode) { // load relevant nodes - auto trkrclusters = findNode::getClass(topNode, "TRKR_CLUSTER"); - if( !trkrclusters ) + auto* trkrclusters = findNode::getClass(topNode, "TRKR_CLUSTER"); + if (!trkrclusters) { std::cout << "MvtxClusterPruner::process_event - TRKR_CLUSTER not found. Doing nothing" << std::endl; return Fun4AllReturnCodes::EVENT_OK; } - auto clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); - if( !clusterhitassoc ) + auto* clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); + if (!clusterhitassoc) { std::cout << "MvtxClusterPruner::process_event - TRKR_CLUSTERHITASSOC not found. Doing nothing" << std::endl; return Fun4AllReturnCodes::EVENT_OK; } // lambda method to create map of cluster keys and associated hits - auto get_cluster_map = [trkrclusters,clusterhitassoc]( TrkrDefs::hitsetkey key ) + auto get_cluster_map = [trkrclusters, clusterhitassoc](TrkrDefs::hitsetkey key) { clustermap_t out; // get all clusters for this hitsetkey - const auto cluster_range= trkrclusters->getClusters(key); - for( const auto& [ckey,cluster]:range_adaptor(cluster_range) ) + const auto cluster_range = trkrclusters->getClusters(key); + for (const auto& [ckey, cluster] : range_adaptor(cluster_range)) { // get associated hits const auto& hit_range = clusterhitassoc->getHits(ckey); hitkeyset_t hitkeys; - std::transform(hit_range.first, hit_range.second, std::inserter(hitkeys,hitkeys.end()), - [](const TrkrClusterHitAssoc::Map::value_type& pair ){ return pair.second; }); - out.emplace(ckey,std::move(hitkeys)); + std::transform(hit_range.first, hit_range.second, std::inserter(hitkeys, hitkeys.end()), + [](const TrkrClusterHitAssoc::Map::value_type& pair) + { return pair.second; }); + out.emplace(ckey, std::move(hitkeys)); } return out; }; // loop over MVTX hitset keys const auto hitsetkeys = trkrclusters->getHitSetKeys(TrkrDefs::mvtxId); - for( const auto& hitsetkey:hitsetkeys ) + for (const auto& hitsetkey : hitsetkeys) { // get layer, stave, chip and current strobe const auto layer = TrkrDefs::getLayer(hitsetkey); @@ -125,111 +133,109 @@ int MvtxClusterPruner::process_event(PHCompositeNode *topNode) const auto cluster_map1 = get_cluster_map(hitsetkey); // get clusters for the next strobe - int next_strobe = current_strobe+1; + int next_strobe = current_strobe + 1; const auto hitsetkey_next_strobe = MvtxDefs::genHitSetKey(layer, stave, chip, next_strobe); const auto clusterk_map2 = get_cluster_map(hitsetkey_next_strobe); // loop over clusters from first range - for( auto [ckey1,hitkeys1]:cluster_map1) + for (auto [ckey1, hitkeys1] : cluster_map1) { // increment counter ++m_cluster_counter_total; // get correcponding cluser - auto cluster1 = Verbosity() ? trkrclusters->findCluster(ckey1):nullptr; + auto* cluster1 = Verbosity() ? trkrclusters->findCluster(ckey1) : nullptr; // loop over clusters from second range - for( auto [ckey2,hitkeys2]:clusterk_map2) + for (auto [ckey2, hitkeys2] : clusterk_map2) { - auto cluster2 = Verbosity() ? trkrclusters->findCluster(ckey2):nullptr; + auto* cluster2 = Verbosity() ? trkrclusters->findCluster(ckey2) : nullptr; - if( m_use_strict_matching ) + if (m_use_strict_matching) { // see if hitsets are identical - if(hitkeys1 == hitkeys2) + if (hitkeys1 == hitkeys2) { // increment counter ++m_cluster_counter_deleted; - if( Verbosity() ) + if (Verbosity()) { std::cout << "Removing cluster "; - print_cluster_information( ckey2, cluster2); + print_cluster_information(ckey2, cluster2); std::cout << "Keeping cluster "; - print_cluster_information( ckey1, cluster1); + print_cluster_information(ckey1, cluster1); } // always remove second cluster trkrclusters->removeCluster(ckey2); break; } - - } else { - + } + else + { // make sure first set is larger than second const bool swapped = hitkeys2.size() > hitkeys1.size(); - if( swapped ) { std::swap(hitkeys2,hitkeys1); } + if (swapped) + { + std::swap(hitkeys2, hitkeys1); + } // see if hitkeys2 is a subset of hitkeys1 - if( std::includes(hitkeys1.begin(), hitkeys1.end(), hitkeys2.begin(), hitkeys2.end()) ) + if (std::includes(hitkeys1.begin(), hitkeys1.end(), hitkeys2.begin(), hitkeys2.end())) { // increment counter ++m_cluster_counter_deleted; - if( swapped ) + if (swapped) { - - if( Verbosity() ) + if (Verbosity()) { std::cout << "Removing cluster "; - print_cluster_information( ckey1, cluster1); + print_cluster_information(ckey1, cluster1); std::cout << "Keeping cluster "; - print_cluster_information( ckey2, cluster2); + print_cluster_information(ckey2, cluster2); } // remove first cluster trkrclusters->removeCluster(ckey1); break; - } else { - - if( Verbosity() ) - { - std::cout << "Removing cluster "; - print_cluster_information( ckey2, cluster2); - - std::cout << "Keeping cluster "; - print_cluster_information( ckey1, cluster1); - } + } + if (Verbosity()) + { + std::cout << "Removing cluster "; + print_cluster_information(ckey2, cluster2); - // remove second cluster - trkrclusters->removeCluster(ckey2); + std::cout << "Keeping cluster "; + print_cluster_information(ckey1, cluster1); } + + // remove second cluster + trkrclusters->removeCluster(ckey2); } - } // strict matching + } // strict matching - } // second cluster loop - } // first cluster loop - } // hitsetkey loop + } // second cluster loop + } // first cluster loop + } // hitsetkey loop return Fun4AllReturnCodes::EVENT_OK; - } //_____________________________________________________________________________ -int MvtxClusterPruner::End(PHCompositeNode * /*topNode*/) +int MvtxClusterPruner::End(PHCompositeNode* /*topNode*/) { - std::cout << "MvtxClusterPruner::End -" - << " m_cluster_counter_total: " << m_cluster_counter_total - << std::endl; - std::cout << "MvtxClusterPruner::End -" - << " m_cluster_counter_deleted: " << m_cluster_counter_deleted - << " fraction: " << double( m_cluster_counter_deleted )/m_cluster_counter_total - << std::endl; + << " m_cluster_counter_total: " << m_cluster_counter_total + << std::endl; + std::cout << "MvtxClusterPruner::End -" + << " m_cluster_counter_deleted: " << m_cluster_counter_deleted + << " fraction: " << double(m_cluster_counter_deleted) / m_cluster_counter_total + << std::endl; return Fun4AllReturnCodes::EVENT_OK; } diff --git a/offline/packages/mvtx/MvtxClusterizer.cc b/offline/packages/mvtx/MvtxClusterizer.cc index 7970b369eb..eb54930acd 100644 --- a/offline/packages/mvtx/MvtxClusterizer.cc +++ b/offline/packages/mvtx/MvtxClusterizer.cc @@ -63,7 +63,7 @@ namespace /// convenience square method template - inline constexpr T square(const T &x) + constexpr T square(const T &x) { return x * x; } @@ -71,65 +71,45 @@ namespace bool MvtxClusterizer::are_adjacent( const std::pair &lhs, - const std::pair &rhs) + const std::pair &rhs) const { if (GetZClustering()) { return - // column adjacent - ( (MvtxDefs::getCol(lhs.first) > MvtxDefs::getCol(rhs.first)) ? - MvtxDefs::getCol(lhs.first)<=MvtxDefs::getCol(rhs.first)+1: - MvtxDefs::getCol(rhs.first)<=MvtxDefs::getCol(lhs.first)+1) && + // column adjacent + ((MvtxDefs::getCol(lhs.first) > MvtxDefs::getCol(rhs.first)) ? MvtxDefs::getCol(lhs.first) <= MvtxDefs::getCol(rhs.first) + 1 : MvtxDefs::getCol(rhs.first) <= MvtxDefs::getCol(lhs.first) + 1) && - // row adjacent - ( (MvtxDefs::getRow(lhs.first) > MvtxDefs::getRow(rhs.first)) ? - MvtxDefs::getRow(lhs.first)<=MvtxDefs::getRow(rhs.first)+1: - MvtxDefs::getRow(rhs.first)<=MvtxDefs::getRow(lhs.first)+1); - - } else { - - return + // row adjacent + ((MvtxDefs::getRow(lhs.first) > MvtxDefs::getRow(rhs.first)) ? MvtxDefs::getRow(lhs.first) <= MvtxDefs::getRow(rhs.first) + 1 : MvtxDefs::getRow(rhs.first) <= MvtxDefs::getRow(lhs.first) + 1); + } + return // column identical - MvtxDefs::getCol(rhs.first)==MvtxDefs::getCol(lhs.first) && + MvtxDefs::getCol(rhs.first) == MvtxDefs::getCol(lhs.first) && // row adjacent - ( (MvtxDefs::getRow(lhs.first) > MvtxDefs::getRow(rhs.first)) ? - MvtxDefs::getRow(lhs.first)<=MvtxDefs::getRow(rhs.first)+1: - MvtxDefs::getRow(rhs.first)<=MvtxDefs::getRow(lhs.first)+1); - - } + ((MvtxDefs::getRow(lhs.first) > MvtxDefs::getRow(rhs.first)) ? MvtxDefs::getRow(lhs.first) <= MvtxDefs::getRow(rhs.first) + 1 : MvtxDefs::getRow(rhs.first) <= MvtxDefs::getRow(lhs.first) + 1); } -bool MvtxClusterizer::are_adjacent(RawHit *lhs, RawHit *rhs) +bool MvtxClusterizer::are_adjacent(RawHit *lhs, RawHit *rhs) const { if (GetZClustering()) { return - // phi adjacent (== column) - ((lhs->getPhiBin() > rhs->getPhiBin()) ? - lhs->getPhiBin() <= rhs->getPhiBin()+1: - rhs->getPhiBin() <= lhs->getPhiBin()+1) && + // phi adjacent (== column) + ((lhs->getPhiBin() > rhs->getPhiBin()) ? lhs->getPhiBin() <= rhs->getPhiBin() + 1 : rhs->getPhiBin() <= lhs->getPhiBin() + 1) && - // time adjacent (== row) - ((lhs->getTBin() > rhs->getTBin()) ? - lhs->getTBin() <= rhs->getTBin()+1: - rhs->getTBin() <= lhs->getTBin()+1); - - } else { - - return + // time adjacent (== row) + ((lhs->getTBin() > rhs->getTBin()) ? lhs->getTBin() <= rhs->getTBin() + 1 : rhs->getTBin() <= lhs->getTBin() + 1); + } + return // phi identical (== column) lhs->getPhiBin() == rhs->getPhiBin() && // time adjacent (== row) - ((lhs->getTBin() > rhs->getTBin()) ? - lhs->getTBin() <= rhs->getTBin()+1: - rhs->getTBin() <= lhs->getTBin()+1); - - } + ((lhs->getTBin() > rhs->getTBin()) ? lhs->getTBin() <= rhs->getTBin() + 1 : rhs->getTBin() <= lhs->getTBin() + 1); } MvtxClusterizer::MvtxClusterizer(const std::string &name) @@ -165,7 +145,7 @@ int MvtxClusterizer::InitRun(PHCompositeNode *topNode) } // Create the Cluster node if required - auto trkrclusters = + auto *trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); if (!trkrclusters) { @@ -183,7 +163,7 @@ int MvtxClusterizer::InitRun(PHCompositeNode *topNode) DetNode->addNode(TrkrClusterContainerNode); } - auto clusterhitassoc = + auto *clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); if (!clusterhitassoc) { @@ -208,14 +188,14 @@ int MvtxClusterizer::InitRun(PHCompositeNode *topNode) if (!mClusHitsVerbose) { PHNodeIterator dstiter(dstNode); - auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); dstNode->addNode(DetNode); } mClusHitsVerbose = new ClusHitsVerbosev1(); - auto newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject"); + auto *newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject"); DetNode->addNode(newNode); } } @@ -227,13 +207,13 @@ int MvtxClusterizer::InitRun(PHCompositeNode *topNode) if (Verbosity() > 0) { std::cout << "====================== MvtxClusterizer::InitRun() " - "=====================" - << std::endl; + "=====================" + << std::endl; std::cout << " Z-dimension Clustering = " << std::boolalpha << m_makeZClustering - << std::noboolalpha << std::endl; + << std::noboolalpha << std::endl; std::cout << "==================================================================" - "=========" - << std::endl; + "=========" + << std::endl; } return Fun4AllReturnCodes::EVENT_OK; @@ -283,7 +263,7 @@ int MvtxClusterizer::process_event(PHCompositeNode *topNode) // reset MVTX clusters and cluster associations const auto hitsetkeys = m_clusterlist->getHitSetKeys(TrkrDefs::mvtxId); - for( const auto& hitsetkey:hitsetkeys) + for (const auto &hitsetkey : hitsetkeys) { m_clusterlist->removeClusters(hitsetkey); m_clusterhitassoc->removeAssocs(hitsetkey); @@ -337,8 +317,8 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) unsigned int chip = MvtxDefs::getChipId(hitsetitr->first); unsigned int strobe = MvtxDefs::getStrobeId(hitsetitr->first); std::cout << "MvtxClusterizer found hitsetkey " << hitsetitr->first - << " layer " << layer << " stave " << stave << " chip " << chip - << " strobe " << strobe << std::endl; + << " layer " << layer << " stave " << stave << " chip " << chip + << " strobe " << strobe << std::endl; } if (Verbosity() > 2) @@ -394,7 +374,7 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) std::vector component(num_vertices(G)); // this is the actual clustering, performed by boost - boost::connected_components(G, &component[0]); + boost::connected_components(G, component.data()); // Loop over the components(hits) compiling a list of the // unique connected groups (ie. clusters). @@ -405,7 +385,7 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) cluster_ids.insert(component[i]); clusters.insert(make_pair(component[i], hitvec[i])); } - for (const auto& clusid:cluster_ids) + for (const auto &clusid : cluster_ids) { auto clusrange = clusters.equal_range(clusid); auto ckey = TrkrDefs::genClusKey(hitset->getHitSetKey(), clusid); @@ -413,7 +393,8 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) // determine the size of the cluster in phi and z std::set phibins; std::set zbins; - std::map m_phi, m_z; // Note, there are no "cut" bins for Svtx Clusters + std::map m_phi; + std::map m_z; // Note, there are no "cut" bins for Svtx Clusters // determine the cluster position... double locxsum = 0.; @@ -426,7 +407,7 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) // we need the geometry object for this layer to get the global positions int layer = TrkrDefs::getLayer(ckey); - auto layergeom = dynamic_cast(geom_container->GetLayerGeom(layer)); + auto *layergeom = dynamic_cast(geom_container->GetLayerGeom(layer)); if (!layergeom) { exit(1); @@ -574,11 +555,11 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) if (Verbosity() > 0) { std::cout << " MvtxClusterizer: cluskey " << ckey << " layer " << layer - << " rad " << layergeom->get_radius() << " phibins " - << phibins.size() << " pitch " << pitch << " phisize " << phisize - << " zbins " << zbins.size() << " length " << length << " zsize " - << zsize << " local x " << locclusx << " local y " << locclusz - << std::endl; + << " rad " << layergeom->get_radius() << " phibins " + << phibins.size() << " pitch " << pitch << " phisize " << phisize + << " zbins " << zbins.size() << " length " << length << " zsize " + << zsize << " local x " << locclusx << " local y " << locclusz + << std::endl; } auto clus = std::make_unique(); @@ -605,7 +586,7 @@ void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) } } // clusitr loop - } // loop over hitsets + } // loop over hitsets if (Verbosity() > 1) { @@ -650,8 +631,8 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) unsigned int chip = MvtxDefs::getChipId(hitsetitr->first); unsigned int strobe = MvtxDefs::getStrobeId(hitsetitr->first); std::cout << "MvtxClusterizer found hitsetkey " << hitsetitr->first - << " layer " << layer << " stave " << stave << " chip " << chip - << " strobe " << strobe << std::endl; + << " layer " << layer << " stave " << stave << " chip " << chip + << " strobe " << strobe << std::endl; } if (Verbosity() > 2) @@ -695,7 +676,7 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) std::vector component(num_vertices(G)); // this is the actual clustering, performed by boost - boost::connected_components(G, &component[0]); + boost::connected_components(G, component.data()); // Loop over the components(hits) compiling a list of the // unique connected groups (ie. clusters). @@ -709,7 +690,7 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) } // std::cout << "found cluster #: "<< clusters.size()<< std::endl; // loop over the componenets and make clusters - for( const auto& clusid:cluster_ids) + for (const auto &clusid : cluster_ids) { auto clusrange = clusters.equal_range(clusid); @@ -731,7 +712,7 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) // we need the geometry object for this layer to get the global positions int layer = TrkrDefs::getLayer(ckey); - auto layergeom = dynamic_cast( + auto *layergeom = dynamic_cast( geom_container->GetLayerGeom(layer)); if (!layergeom) { @@ -845,11 +826,11 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) if (Verbosity() > 0) { std::cout << " MvtxClusterizer: cluskey " << ckey << " layer " << layer - << " rad " << layergeom->get_radius() << " phibins " - << phibins.size() << " pitch " << pitch << " phisize " << phisize - << " zbins " << zbins.size() << " length " << length << " zsize " - << zsize << " local x " << locclusx << " local y " << locclusz - << std::endl; + << " rad " << layergeom->get_radius() << " phibins " + << phibins.size() << " pitch " << pitch << " phisize " << phisize + << " zbins " << zbins.size() << " length " << length << " zsize " + << zsize << " local x " << locclusx << " local y " << locclusz + << std::endl; } auto clus = std::make_unique(); @@ -875,7 +856,7 @@ void MvtxClusterizer::ClusterMvtxRaw(PHCompositeNode *topNode) m_clusterlist->addClusterSpecifyKey(ckey, clus.release()); } } // clusitr loop - } // loop over hitsets + } // loop over hitsets if (Verbosity() > 1) { @@ -898,11 +879,11 @@ void MvtxClusterizer::PrintClusters(PHCompositeNode *topNode) } std::cout << "================= After MvtxClusterizer::process_event() " - "====================" - << std::endl; + "====================" + << std::endl; std::cout << " There are " << clusterlist->size() - << " clusters recorded: " << std::endl; + << " clusters recorded: " << std::endl; if (Verbosity() > 3) { @@ -910,8 +891,8 @@ void MvtxClusterizer::PrintClusters(PHCompositeNode *topNode) } std::cout << "==================================================================" - "=========" - << std::endl; + "=========" + << std::endl; } return; diff --git a/offline/packages/mvtx/MvtxClusterizer.h b/offline/packages/mvtx/MvtxClusterizer.h index ce4794b5d1..10c65914cd 100644 --- a/offline/packages/mvtx/MvtxClusterizer.h +++ b/offline/packages/mvtx/MvtxClusterizer.h @@ -65,8 +65,8 @@ class MvtxClusterizer : public SubsysReco private: // bool are_adjacent(const pixel lhs, const pixel rhs); bool record_ClusHitsVerbose{false}; - bool are_adjacent(const std::pair &lhs, const std::pair &rhs); - bool are_adjacent(RawHit *lhs, RawHit *rhs); + bool are_adjacent(const std::pair &lhs, const std::pair &rhs) const; + bool are_adjacent(RawHit *lhs, RawHit *rhs) const; void ClusterMvtx(PHCompositeNode *topNode); void ClusterMvtxRaw(PHCompositeNode *topNode); diff --git a/offline/packages/mvtx/MvtxHitPruner.cc b/offline/packages/mvtx/MvtxHitPruner.cc index 4f141a4c3f..58cbb423c4 100644 --- a/offline/packages/mvtx/MvtxHitPruner.cc +++ b/offline/packages/mvtx/MvtxHitPruner.cc @@ -51,26 +51,33 @@ namespace { //! range adaptor to be able to use range-based for loop - template class range_adaptor + template + class range_adaptor { - public: - range_adaptor( const T& range ):m_range(range){} - const typename T::first_type& begin() {return m_range.first;} - const typename T::second_type& end() {return m_range.second;} - private: + public: + explicit range_adaptor(const T& range) + : m_range(range) + { + } + const typename T::first_type& begin() { return m_range.first; } + const typename T::second_type& end() { return m_range.second; } + + private: T m_range; }; -} +} // namespace -MvtxHitPruner::MvtxHitPruner(const std::string &name) +MvtxHitPruner::MvtxHitPruner(const std::string& name) : SubsysReco(name) { } -int MvtxHitPruner::InitRun(PHCompositeNode * /*topNode*/) -{ return Fun4AllReturnCodes::EVENT_OK; } +int MvtxHitPruner::InitRun(PHCompositeNode* /*topNode*/) +{ + return Fun4AllReturnCodes::EVENT_OK; +} -int MvtxHitPruner::process_event(PHCompositeNode *topNode) +int MvtxHitPruner::process_event(PHCompositeNode* topNode) { // get node containing the digitized hits m_hits = findNode::getClass(topNode, "TRKR_HITSET"); @@ -93,12 +100,14 @@ int MvtxHitPruner::process_event(PHCompositeNode *topNode) std::set bare_hitset_set; const auto hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId); - for( const auto& [hitsetkey,hitset]:range_adaptor(hitsetrange) ) + for (const auto& [hitsetkey, hitset] : range_adaptor(hitsetrange)) { - // get strobe, skip if already zero const int strobe = MvtxDefs::getStrobeId(hitsetkey); - if( strobe == 0 ) continue; + if (strobe == 0) + { + continue; + } // get the hitsetkey value for strobe 0 const auto bare_hitsetkey = MvtxDefs::resetStrobe(hitsetkey); @@ -117,43 +126,46 @@ int MvtxHitPruner::process_event(PHCompositeNode *topNode) for (const auto& bare_hitsetkey : bare_hitset_set) { // find matching hitset of creater - auto bare_hitset = (m_hits->findOrAddHitSet(bare_hitsetkey))->second; + auto* bare_hitset = (m_hits->findOrAddHitSet(bare_hitsetkey))->second; if (Verbosity()) { std::cout - << "MvtxHitPruner::process_event - bare_hitset " << bare_hitsetkey - << " initially has " << bare_hitset->size() << " hits " - << std::endl; + << "MvtxHitPruner::process_event - bare_hitset " << bare_hitsetkey + << " initially has " << bare_hitset->size() << " hits " + << std::endl; } // get all hitsets with non-zero strobe that match the bare hitset key auto bare_hitsetrange = hitset_multimap.equal_range(bare_hitsetkey); - for( const auto& [unused,hitsetkey]:range_adaptor(bare_hitsetrange) ) + for (const auto& [unused, hitsetkey] : range_adaptor(bare_hitsetrange)) { const int strobe = MvtxDefs::getStrobeId(hitsetkey); - if( strobe == 0 ) continue; + if (strobe == 0) + { + continue; + } if (Verbosity()) { std::cout << "MvtxHitPruner::process_event -" - << " process hitsetkey " << hitsetkey - << " from strobe " << strobe - << " for bare_hitsetkey " << bare_hitsetkey - << std::endl; + << " process hitsetkey " << hitsetkey + << " from strobe " << strobe + << " for bare_hitsetkey " << bare_hitsetkey + << std::endl; } // copy all hits to the hitset with strobe 0 - auto hitset = m_hits->findHitSet(hitsetkey); + auto* hitset = m_hits->findHitSet(hitsetkey); if (Verbosity()) { std::cout << "MvtxHitPruner::process_event - hitsetkey " << hitsetkey - << " has strobe " << strobe << " and has " << hitset->size() - << " hits, so copy it" << std::endl; + << " has strobe " << strobe << " and has " << hitset->size() + << " hits, so copy it" << std::endl; } TrkrHitSet::ConstRange hitrangei = hitset->getHits(); - for( const auto& [hitkey,old_hit]:range_adaptor(hitrangei) ) + for (const auto& [hitkey, old_hit] : range_adaptor(hitrangei)) { if (Verbosity()) { @@ -166,9 +178,9 @@ int MvtxHitPruner::process_event(PHCompositeNode *topNode) if (Verbosity()) { std::cout - << "MvtxHitPruner::process_event - hitkey " << hitkey - << " is already in bare hitsest, do not copy" - << std::endl; + << "MvtxHitPruner::process_event - hitkey " << hitkey + << " is already in bare hitsest, do not copy" + << std::endl; } continue; } @@ -177,11 +189,11 @@ int MvtxHitPruner::process_event(PHCompositeNode *topNode) if (Verbosity()) { std::cout - << "MvtxHitPruner::process_event - copying over hitkey " - << hitkey << std::endl; + << "MvtxHitPruner::process_event - copying over hitkey " + << hitkey << std::endl; } - auto new_hit = new TrkrHitv2; + auto* new_hit = new TrkrHitv2; new_hit->CopyFrom(old_hit); bare_hitset->addHitSpecificKey(hitkey, new_hit); } diff --git a/offline/packages/mvtx/SegmentationAlpide.cc b/offline/packages/mvtx/SegmentationAlpide.cc index 6b9aecddd0..0b1aae48cf 100644 --- a/offline/packages/mvtx/SegmentationAlpide.cc +++ b/offline/packages/mvtx/SegmentationAlpide.cc @@ -5,16 +5,15 @@ */ #include "SegmentationAlpide.h" -#include - #include +#include void SegmentationAlpide::print() { - std::cout << (boost::format("Pixel size: %.2f (along %d rows) %.2f (along %d columns) microns") % (PitchRow * 1e4) % NRows % (PitchCol * 1e4) % NCols).str() + std::cout << std::format("Pixel size: {:.2f} (along {} rows) {:.2f} (along {} columns) microns", (PitchRow * 1e4), NRows, (PitchCol * 1e4), NCols) << std::endl; - std::cout << (boost::format("Passive edges: bottom: %.2f, top: %.2f, left/right: %.2f microns") % (PassiveEdgeReadOut * 1e4) % (PassiveEdgeTop * 1e4) % (PassiveEdgeSide * 1e4)).str() + std::cout << std::format("Passive edges: bottom: {:.2f}, top: {:.2f}, left/right: {:.2f} microns", (PassiveEdgeReadOut * 1e4), (PassiveEdgeTop * 1e4), (PassiveEdgeSide * 1e4)) << std::endl; - std::cout << (boost::format("Active/Total size: %.6f/%.6f (rows) %.6f/%.6f (cols) cm") % ActiveMatrixSizeRows % SensorSizeRows % ActiveMatrixSizeCols % SensorSizeCols).str() + std::cout << std::format("Active/Total size: {:.6f}/{:.6f} (rows) {:.6f}/{:.6f} (cols) cm", ActiveMatrixSizeRows, SensorSizeRows, ActiveMatrixSizeCols, SensorSizeCols) << std::endl; } diff --git a/offline/packages/tpc/LaserClusterizer.cc b/offline/packages/tpc/LaserClusterizer.cc index 170f95944d..754ee8d248 100644 --- a/offline/packages/tpc/LaserClusterizer.cc +++ b/offline/packages/tpc/LaserClusterizer.cc @@ -31,13 +31,10 @@ #include #include +#include #include #include #include -//#include -//#include -#include -//#include #include #include @@ -46,15 +43,16 @@ #include #include // for sqrt, cos, sin +#include #include #include #include // for _Rb_tree_cons... +#include +#include #include +#include #include // for pair #include -#include -#include -#include #include @@ -68,12 +66,9 @@ using adcKey = std::pair; using pointKeyLaser = std::pair; using hitData = std::pair; - -int layerMins[3] = {7,23,39}; +int layerMins[3] = {7, 23, 39}; int layerMaxes[3] = {22, 38, 54}; - - namespace { struct thread_data @@ -101,28 +96,25 @@ namespace pthread_mutex_t mythreadlock; const std::vector neighborOffsets = { - point(1, 0, 0), point(-1, 0, 0), - point(0, 1, 0), point(0, -1, 0), - point(0, 0, 1), point(0, 0, -1), - point(0, 0, 2), point(0, 0, -2) - }; - + point(1, 0, 0), point(-1, 0, 0), + point(0, 1, 0), point(0, -1, 0), + point(0, 0, 1), point(0, 0, -1), + point(0, 0, 2), point(0, 0, -2)}; - double layerFunction(double *x, double *par) + double layerFunction(double *x, const double *par) { double A = par[0]; double mu = par[1]; - double binCenter = round(x[0]); double overlapLow = std::max(binCenter - 0.5, mu - 0.5); double overlapHigh = std::min(binCenter + 0.5, mu + 0.5); double overlap = overlapHigh - overlapLow; - if(overlap <= 0.0) + if (overlap <= 0.0) { return 0.0; } - return A*overlap; + return A * overlap; /* if(fabs(x[0] - mu) < 1) { @@ -131,28 +123,26 @@ namespace } return 0.0; */ - - } double phiFunction(double *x, double *par) { - if(par[2] < 0.0) + if (par[2] < 0.0) { return 0.0; } - return par[0] * TMath::Gaus(x[0],par[1],par[2],false); + return par[0] * TMath::Gaus(x[0], par[1], par[2], false); } double timeFunction(double *x, double *par) { - if(par[2] < 0.0) + if (par[2] < 0.0) { return 0.0; } - double g = TMath::Gaus(x[0],par[1],par[2],true); - double cdf = 1 + TMath::Erfc(par[3]*(x[0]-par[1])/(sqrt(2.0)*par[2])); - return par[0]*g*cdf; + double g = TMath::Gaus(x[0], par[1], par[2], true); + double cdf = 1 + TMath::Erfc(par[3] * (x[0] - par[1]) / (sqrt(2.0) * par[2])); + return par[0] * g * cdf; } void findConnectedRegions3(std::vector &clusHits, std::pair &maxKey) @@ -160,21 +150,22 @@ namespace std::vector> regions; std::vector unvisited; - for(auto &clusHit : clusHits) + unvisited.reserve(clusHits.size()); + for (auto &clusHit : clusHits) { unvisited.push_back(clusHit); } - while(!unvisited.empty()) + while (!unvisited.empty()) { std::vector region; std::queue q; unsigned int mIndex = 0; - int i=0; - for(auto hit : unvisited) + int i = 0; + for (auto hit : unvisited) { - if(hit.second.second.first == maxKey.first && hit.second.second.second == maxKey.second) + if (hit.second.second.first == maxKey.first && hit.second.second.second == maxKey.second) { mIndex = i; break; @@ -183,11 +174,11 @@ namespace } auto seed = unvisited[mIndex]; - unvisited.erase(unvisited.begin()+mIndex); + unvisited.erase(unvisited.begin() + mIndex); q.push(seed); region.push_back(seed); - while(!q.empty()) + while (!q.empty()) { float ix = q.front().first.get<0>(); float iy = q.front().first.get<1>(); @@ -200,13 +191,12 @@ namespace float ny = iy + neigh.get<1>(); float nz = iz + neigh.get<2>(); - for(unsigned int v=0; v() - nx) < 0.01 && fabs(unvisited[v].first.get<1>() - ny) < 0.01 && fabs(unvisited[v].first.get<2>() - nz) < 0.01) + if (fabs(unvisited[v].first.get<0>() - nx) < 0.01 && fabs(unvisited[v].first.get<1>() - ny) < 0.01 && fabs(unvisited[v].first.get<2>() - nz) < 0.01) { auto newSeed = unvisited[v]; - unvisited.erase(unvisited.begin()+v); + unvisited.erase(unvisited.begin() + v); q.push(newSeed); region.push_back(newSeed); break; @@ -215,18 +205,15 @@ namespace } } regions.push_back(region); - } clusHits.clear(); - for(auto hit : regions[0]) + for (auto hit : regions[0]) { clusHits.push_back(hit); } - } - void remove_hits(std::vector &clusHits, bgi::rtree> &rtree, std::multimap &adcMap) { for (auto &clusHit : clusHits) @@ -237,51 +224,44 @@ namespace for (auto iterAdc = adcMap.begin(); iterAdc != adcMap.end();) { - if(iterAdc->second.second == spechitkey) - { - iterAdc = adcMap.erase(iterAdc); - break; - } - else - { - ++iterAdc; - } + if (iterAdc->second.second == spechitkey) + { + iterAdc = adcMap.erase(iterAdc); + break; + } + + ++iterAdc; } } - } void calc_cluster_parameter(std::vector &clusHits, thread_data &my_data, std::pair maxADCKey) { - - - findConnectedRegions3(clusHits, maxADCKey); - double rSum = 0.0; double phiSum = 0.0; double tSum = 0.0; - + double layerSum = 0.0; double iphiSum = 0.0; double itSum = 0.0; - + double adcSum = 0.0; - + double maxAdc = 0.0; TrkrDefs::hitsetkey maxKey = 0; - + unsigned int nHits = clusHits.size(); - + auto *clus = new LaserClusterv2; - + int meanSide = 0; - + std::vector usedLayer; std::vector usedIPhi; std::vector usedIT; - + double meanLayer = 0.0; double meanIPhi = 0.0; double meanIT = 0.0; @@ -293,113 +273,110 @@ namespace unsigned int adc = clusHit.second.first; int side = TpcDefs::getSide(spechitkey.second); - + if (side) { - meanSide++; + meanSide++; } else { - meanSide--; + meanSide--; } - + PHG4TpcGeom *layergeom = my_data.geom_container->GetLayerCellGeom((int) coords[0]); - + double r = layergeom->get_radius(); double phi = layergeom->get_phi(coords[1], side); double t = layergeom->get_zcenter(fabs(coords[2])); - + double hitzdriftlength = t * my_data.tGeometry->get_drift_velocity(); double hitZ = my_data.tdriftmax * my_data.tGeometry->get_drift_velocity() - hitzdriftlength; - - - bool foundLayer = false; - for (float i : usedLayer) - { - if (coords[0] == i) - { - foundLayer = true; - break; - } - } - - if (!foundLayer) - { - usedLayer.push_back(coords[0]); - } - - bool foundIPhi = false; - for (float i : usedIPhi) - { - if (coords[1] == i) - { - foundIPhi = true; - break; - } - } - - if (!foundIPhi) - { - usedIPhi.push_back(coords[1]); - } - - bool foundIT = false; - for (float i : usedIT) - { - if (coords[2] == i) - { - foundIT = true; - break; - } - } - - if (!foundIT) - { - usedIT.push_back(coords[2]); - } - - clus->addHit(); - clus->setHitLayer(clus->getNhits() - 1, coords[0]); - clus->setHitIPhi(clus->getNhits() - 1, coords[1]); - clus->setHitIT(clus->getNhits() - 1, coords[2]); - clus->setHitX(clus->getNhits() - 1, r * cos(phi)); - clus->setHitY(clus->getNhits() - 1, r * sin(phi)); - clus->setHitZ(clus->getNhits() - 1, hitZ); - clus->setHitAdc(clus->getNhits() - 1, (float) adc); - - rSum += r * adc; - phiSum += phi * adc; - tSum += t * adc; - - layerSum += coords[0] * adc; - iphiSum += coords[1] * adc; - itSum += coords[2] * adc; - - meanLayer += coords[0]; - meanIPhi += coords[1]; - meanIT += coords[2]; - - adcSum += adc; - - if (adc > maxAdc) - { - maxAdc = adc; - maxKey = spechitkey.second; - } - - } - + + bool foundLayer = false; + for (float i : usedLayer) + { + if (coords[0] == i) + { + foundLayer = true; + break; + } + } + + if (!foundLayer) + { + usedLayer.push_back(coords[0]); + } + + bool foundIPhi = false; + for (float i : usedIPhi) + { + if (coords[1] == i) + { + foundIPhi = true; + break; + } + } + + if (!foundIPhi) + { + usedIPhi.push_back(coords[1]); + } + + bool foundIT = false; + for (float i : usedIT) + { + if (coords[2] == i) + { + foundIT = true; + break; + } + } + + if (!foundIT) + { + usedIT.push_back(coords[2]); + } + + clus->addHit(); + clus->setHitLayer(clus->getNhits() - 1, coords[0]); + clus->setHitIPhi(clus->getNhits() - 1, coords[1]); + clus->setHitIT(clus->getNhits() - 1, coords[2]); + clus->setHitX(clus->getNhits() - 1, r * cos(phi)); + clus->setHitY(clus->getNhits() - 1, r * sin(phi)); + clus->setHitZ(clus->getNhits() - 1, hitZ); + clus->setHitAdc(clus->getNhits() - 1, (float) adc); + + rSum += r * adc; + phiSum += phi * adc; + tSum += t * adc; + + layerSum += coords[0] * adc; + iphiSum += coords[1] * adc; + itSum += coords[2] * adc; + + meanLayer += coords[0]; + meanIPhi += coords[1]; + meanIT += coords[2]; + + adcSum += adc; + + if (adc > maxAdc) + { + maxAdc = adc; + maxKey = spechitkey.second; + } + } + if (nHits == 0) { return; } - double clusR = rSum / adcSum; double clusPhi = phiSum / adcSum; double clusT = tSum / adcSum; double zdriftlength = clusT * my_data.tGeometry->get_drift_velocity(); - + double clusX = clusR * cos(clusPhi); double clusY = clusR * sin(clusPhi); double clusZ = my_data.tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength; @@ -408,10 +385,10 @@ namespace clusZ = -clusZ; for (int i = 0; i < (int) clus->getNhits(); i++) { - clus->setHitZ(i, -1 * clus->getHitZ(i)); + clus->setHitZ(i, -1 * clus->getHitZ(i)); } } - + std::sort(usedLayer.begin(), usedLayer.end()); std::sort(usedIPhi.begin(), usedIPhi.end()); std::sort(usedIT.begin(), usedIT.end()); @@ -419,30 +396,28 @@ namespace meanLayer = meanLayer / nHits; meanIPhi = meanIPhi / nHits; meanIT = meanIT / nHits; - + double sigmaLayer = 0.0; double sigmaIPhi = 0.0; double sigmaIT = 0.0; - + double sigmaWeightedLayer = 0.0; double sigmaWeightedIPhi = 0.0; double sigmaWeightedIT = 0.0; - - pthread_mutex_lock(&mythreadlock); - my_data.hitHist = new TH3D(Form("hitHist_event%d_side%d_sector%d_module%d_cluster%d",my_data.eventNum,(int)my_data.side,(int)my_data.sector,(int)my_data.module,(int)my_data.cluster_vector.size()),";layer;iphi;it",usedLayer.size()+2,usedLayer[0]-1.5,*usedLayer.rbegin()+1.5,usedIPhi.size()+2,usedIPhi[0]-1.5,*usedIPhi.rbegin()+1.5,usedIT.size()+2,usedIT[0]-1.5,*usedIT.rbegin()+1.5); + pthread_mutex_lock(&mythreadlock); + my_data.hitHist = new TH3D(std::format("hitHist_event{}_side{}_sector{}_module{}_cluster{}", my_data.eventNum, (int) my_data.side, (int) my_data.sector, (int) my_data.module, (int) my_data.cluster_vector.size()).c_str(), ";layer;iphi;it", usedLayer.size() + 2, usedLayer[0] - 1.5, *usedLayer.rbegin() + 1.5, usedIPhi.size() + 2, usedIPhi[0] - 1.5, *usedIPhi.rbegin() + 1.5, usedIT.size() + 2, usedIT[0] - 1.5, *usedIT.rbegin() + 1.5); - //TH3D *hitHist = new TH3D(Form("hitHist_event%d_side%d_sector%d_module%d_cluster%d",my_data.eventNum,(int)my_data.side,(int)my_data.sector,(int)my_data.module,(int)my_data.cluster_vector.size()),";layer;iphi;it",usedLayer.size()+2,usedLayer[0]-1.5,*usedLayer.rbegin()+1.5,usedIPhi.size()+2,usedIPhi[0]-1.5,*usedIPhi.rbegin()+1.5,usedIT.size()+2,usedIT[0]-1.5,*usedIT.rbegin()+1.5); + // TH3D *hitHist = new TH3D(Form("hitHist_event%d_side%d_sector%d_module%d_cluster%d",my_data.eventNum,(int)my_data.side,(int)my_data.sector,(int)my_data.module,(int)my_data.cluster_vector.size()),";layer;iphi;it",usedLayer.size()+2,usedLayer[0]-1.5,*usedLayer.rbegin()+1.5,usedIPhi.size()+2,usedIPhi[0]-1.5,*usedIPhi.rbegin()+1.5,usedIT.size()+2,usedIT[0]-1.5,*usedIT.rbegin()+1.5); for (int i = 0; i < (int) clus->getNhits(); i++) { - my_data.hitHist->Fill(clus->getHitLayer(i), clus->getHitIPhi(i), clus->getHitIT(i), clus->getHitAdc(i)); sigmaLayer += pow(clus->getHitLayer(i) - meanLayer, 2); sigmaIPhi += pow(clus->getHitIPhi(i) - meanIPhi, 2); sigmaIT += pow(clus->getHitIT(i) - meanIT, 2); - + sigmaWeightedLayer += clus->getHitAdc(i) * pow(clus->getHitLayer(i) - (layerSum / adcSum), 2); sigmaWeightedIPhi += clus->getHitAdc(i) * pow(clus->getHitIPhi(i) - (iphiSum / adcSum), 2); sigmaWeightedIT += clus->getHitAdc(i) * pow(clus->getHitIT(i) - (itSum / adcSum), 2); @@ -451,20 +426,18 @@ namespace bool fitSuccess = false; ROOT::Fit::Fitter *fit3D = new ROOT::Fit::Fitter; - if(my_data.doFitting) + if (my_data.doFitting) { - double par_init[7] = { - maxAdc, - meanLayer, - meanIPhi, 0.75, - meanIT, 0.5, 1 - }; + maxAdc, + meanLayer, + meanIPhi, 0.75, + meanIT, 0.5, 1}; double satThreshold = 900.0; double sigma_ADC = 20.0; - auto nll = [&](const double* par) + auto nll = [&](const double *par) { double nll_val = 0.0; @@ -472,37 +445,37 @@ namespace int ny = my_data.hitHist->GetNbinsY(); int nz = my_data.hitHist->GetNbinsZ(); - double parLayer[2] = {1.0,par[1]}; - double parPhi[4] = {1.0,par[2],par[3]}; - double parTime[4] = {1.0,par[4],par[5],par[6]}; + double parLayer[2] = {1.0, par[1]}; + double parPhi[4] = {1.0, par[2], par[3]}; + double parTime[4] = {1.0, par[4], par[5], par[6]}; double xyz[3]; for (int i = 1; i <= nx; ++i) { - xyz[0] = my_data.hitHist->GetXaxis()->GetBinCenter(i); + xyz[0] = my_data.hitHist->GetXaxis()->GetBinCenter(i); for (int j = 1; j <= ny; ++j) { - xyz[1] = my_data.hitHist->GetYaxis()->GetBinCenter(j); + xyz[1] = my_data.hitHist->GetYaxis()->GetBinCenter(j); for (int k = 1; k <= nz; ++k) { xyz[2] = my_data.hitHist->GetZaxis()->GetBinCenter(k); double observed = my_data.hitHist->GetBinContent(i, j, k); - double expected = par[0]*layerFunction(&xyz[0], parLayer)*phiFunction(&xyz[1], parPhi)*timeFunction(&xyz[2], parTime); + double expected = par[0] * layerFunction(&xyz[0], parLayer) * phiFunction(&xyz[1], parPhi) * timeFunction(&xyz[2], parTime); - if(observed <= my_data.adc_threshold) + if (observed <= my_data.adc_threshold) { double arg = (expected - my_data.adc_threshold) / (sqrt(2.0) * sigma_ADC); double tail_prob = 0.5 * TMath::Erfc(arg); nll_val -= log(tail_prob + 1e-12); } - else if(observed < satThreshold) + else if (observed < satThreshold) { double resid = (observed - expected) / sigma_ADC; nll_val += 0.5 * (resid * resid + log(2 * TMath::Pi() * sigma_ADC * sigma_ADC)); } - else if(observed >= satThreshold) + else if (observed >= satThreshold) { double arg = (satThreshold - expected) / (sqrt(2.0) * sigma_ADC); double tail_prob = 0.5 * TMath::Erfc(arg); @@ -518,35 +491,33 @@ namespace fit3D->Config().ParSettings(0).SetName("amp"); fit3D->Config().ParSettings(0).SetStepSize(10); - fit3D->Config().ParSettings(0).SetLimits(0,5000); + fit3D->Config().ParSettings(0).SetLimits(0, 5000); fit3D->Config().ParSettings(1).SetName("mu_layer"); fit3D->Config().ParSettings(1).SetStepSize(0.1); - fit3D->Config().ParSettings(1).SetLimits(usedLayer[0],*usedLayer.rbegin()); + fit3D->Config().ParSettings(1).SetLimits(usedLayer[0], *usedLayer.rbegin()); fit3D->Config().ParSettings(2).SetName("mu_phi"); fit3D->Config().ParSettings(2).SetStepSize(0.1); - fit3D->Config().ParSettings(2).SetLimits(usedIPhi[0],*usedIPhi.rbegin()); + fit3D->Config().ParSettings(2).SetLimits(usedIPhi[0], *usedIPhi.rbegin()); fit3D->Config().ParSettings(3).SetName("sig_phi"); fit3D->Config().ParSettings(3).SetStepSize(0.1); - fit3D->Config().ParSettings(3).SetLimits(0.01,2); + fit3D->Config().ParSettings(3).SetLimits(0.01, 2); fit3D->Config().ParSettings(4).SetName("mu_t"); fit3D->Config().ParSettings(4).SetStepSize(0.1); - fit3D->Config().ParSettings(4).SetLimits(usedIT[0],*usedIT.rbegin()); + fit3D->Config().ParSettings(4).SetLimits(usedIT[0], *usedIT.rbegin()); fit3D->Config().ParSettings(5).SetName("sig_t"); fit3D->Config().ParSettings(5).SetStepSize(0.1); - fit3D->Config().ParSettings(5).SetLimits(0.01,10); + fit3D->Config().ParSettings(5).SetLimits(0.01, 10); fit3D->Config().ParSettings(6).SetName("lambda_t"); fit3D->Config().ParSettings(6).SetStepSize(0.01); - fit3D->Config().ParSettings(6).SetLimits(0,5); + fit3D->Config().ParSettings(6).SetLimits(0, 5); - - if(usedLayer.size() == 1) + if (usedLayer.size() == 1) { fit3D->Config().ParSettings(1).Fix(); } fitSuccess = fit3D->FitFCN(); - if (my_data.Verbosity > 2) { std::cout << "fit success: " << fitSuccess << std::endl; @@ -554,13 +525,9 @@ namespace } pthread_mutex_unlock(&mythreadlock); - - - if(my_data.doFitting && fitSuccess) + if (my_data.doFitting && fitSuccess) { - - const ROOT::Fit::FitResult& result = fit3D->Result(); - + const ROOT::Fit::FitResult &result = fit3D->Result(); PHG4TpcGeom *layergeomLow = my_data.geom_container->GetLayerCellGeom((int) floor(result.Parameter(1))); PHG4TpcGeom *layergeomHigh = my_data.geom_container->GetLayerCellGeom((int) ceil(result.Parameter(1))); @@ -569,12 +536,12 @@ namespace double RHigh = layergeomHigh->get_radius(); double phiHigh_RLow = -999.0; - if(ceil(result.Parameter(2)) < layergeomLow->get_phibins()) + if (ceil(result.Parameter(2)) < layergeomLow->get_phibins()) { phiHigh_RLow = layergeomLow->get_phi(ceil(result.Parameter(2)), (meanSide < 0 ? 0 : 1)); } double phiHigh_RHigh = -999.0; - if(ceil(result.Parameter(2)) < layergeomHigh->get_phibins()) + if (ceil(result.Parameter(2)) < layergeomHigh->get_phibins()) { phiHigh_RHigh = layergeomHigh->get_phi(ceil(result.Parameter(2)), (meanSide < 0 ? 0 : 1)); } @@ -587,18 +554,17 @@ namespace double meanPhi_RLow = ((result.Parameter(2) - floor(result.Parameter(2)))) * (phiHigh_RLow - phiLow_RLow) + phiLow_RLow; double meanPhi_RHigh = ((result.Parameter(2) - floor(result.Parameter(2)))) * (phiHigh_RHigh - phiLow_RHigh) + phiLow_RHigh; - double meanPhi = 0.5*(meanPhi_RLow + meanPhi_RHigh); - if(phiHigh_RLow == -999.0 && phiHigh_RHigh != -999.0) + double meanPhi = 0.5 * (meanPhi_RLow + meanPhi_RHigh); + if (phiHigh_RLow == -999.0 && phiHigh_RHigh != -999.0) { meanPhi = meanPhi_RHigh; } - else if(phiHigh_RLow != -999.0 && phiHigh_RHigh == -999.0) + else if (phiHigh_RLow != -999.0 && phiHigh_RHigh == -999.0) { meanPhi = meanPhi_RLow; } - - if(phiHigh_RLow == -999.0 && phiHigh_RHigh == -999.0) + if (phiHigh_RLow == -999.0 && phiHigh_RHigh == -999.0) { clus->setAdc(adcSum); clus->setX(clusX); @@ -619,10 +585,10 @@ namespace clus->setSDWeightedIT(sqrt(sigmaWeightedIT / adcSum)); } else - { + { clus->setAdc(adcSum); - clus->setX(meanR*cos(meanPhi)); - clus->setY(meanR*sin(meanPhi)); + clus->setX(meanR * cos(meanPhi)); + clus->setY(meanR * sin(meanPhi)); clus->setZ(clusZ); clus->setFitMode(true); clus->setLayer(result.Parameter(1)); @@ -663,25 +629,18 @@ namespace const auto ckey = TrkrDefs::genClusKey(maxKey, my_data.cluster_vector.size()); my_data.cluster_vector.push_back(clus); my_data.cluster_key_vector.push_back(ckey); - - if(fit3D) - { - delete fit3D; - } - if(my_data.hitHist) + delete fit3D; + + if (my_data.hitHist) { delete my_data.hitHist; my_data.hitHist = nullptr; } - - } - void ProcessModuleData(thread_data *my_data) { - if (my_data->Verbosity > 2) { pthread_mutex_lock(&mythreadlock); @@ -693,85 +652,84 @@ namespace std::multimap adcMap; - if (my_data->hitsets.size() == 0) + if (my_data->hitsets.empty()) { return; } - for(int i=0; i<(int)my_data->hitsets.size(); i++) + for (int i = 0; i < (int) my_data->hitsets.size(); i++) { auto *hitset = my_data->hitsets[i]; unsigned int layer = my_data->layers[i]; bool side = my_data->side; unsigned int sector = my_data->sector; - TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, (int)side); + TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, (int) side); TrkrHitSet::ConstRange hitrangei = hitset->getHits(); for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr) { - float_t fadc = hitr->second->getAdc(); - unsigned short adc = 0; - if (fadc > my_data->adc_threshold) - { - adc = (unsigned short) fadc; - } - else - { - continue; - } - - int iphi = TpcDefs::getPad(hitr->first); - int it = TpcDefs::getTBin(hitr->first); - - if(fabs(it - my_data->peakTimeBin) > 5) - { - continue; - } - - point coords = point((int) layer, iphi, it); - - std::vector testduplicate; - rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001), - point(layer + 0.001, iphi + 0.001, it + 0.001))), - std::back_inserter(testduplicate)); - if (!testduplicate.empty()) - { - testduplicate.clear(); - continue; - } - - TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it); - - auto spechitkey = std::make_pair(hitKey, hitsetKey); - pointKeyLaser coordsKey = std::make_pair(coords, spechitkey); - adcMap.insert(std::make_pair(adc, coordsKey)); + float_t fadc = hitr->second->getAdc(); + unsigned short adc = 0; + if (fadc > my_data->adc_threshold) + { + adc = (unsigned short) fadc; + } + else + { + continue; + } + + int iphi = TpcDefs::getPad(hitr->first); + int it = TpcDefs::getTBin(hitr->first); + + if (fabs(it - my_data->peakTimeBin) > 5) + { + continue; + } + + point coords = point((int) layer, iphi, it); + + std::vector testduplicate; + rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001), + point(layer + 0.001, iphi + 0.001, it + 0.001))), + std::back_inserter(testduplicate)); + if (!testduplicate.empty()) + { + testduplicate.clear(); + continue; + } + + TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it); + + auto spechitkey = std::make_pair(hitKey, hitsetKey); + pointKeyLaser coordsKey = std::make_pair(coords, spechitkey); + adcMap.insert(std::make_pair(adc, coordsKey)); auto adckey = std::make_pair(adc, spechitkey); - rtree.insert(std::make_pair(point(1.0*layer, 1.0*iphi, 1.0*it), adckey)); + rtree.insert(std::make_pair(point(1.0 * layer, 1.0 * iphi, 1.0 * it), adckey)); } } - //finished filling rtree + // finished filling rtree - while (adcMap.size() > 0) + while (!adcMap.empty()) { auto iterKey = adcMap.rbegin(); - if(iterKey == adcMap.rend()) + if (iterKey == adcMap.rend()) { - break; + break; } - auto coords = iterKey->second.first; int layer = coords.get<0>(); int iphi = coords.get<1>(); int it = coords.get<2>(); - + if (my_data->Verbosity > 2) { pthread_mutex_lock(&mythreadlock); - std::cout << "working on cluster " << my_data->cluster_vector.size() << " side: " << my_data->side << " sector: " << my_data->sector << " module: " << (layer<23 ? 1 : (layer<39 ? 2 : 3) ) << std::endl; + // NOLINTNEXTLINE (readability-avoid-nested-conditional-operator) + std::cout << "working on cluster " << my_data->cluster_vector.size() << " side: " << my_data->side << " sector: " << my_data->sector << " module: " << (layer < 23 ? 1 : (layer < 39 ? 2 : 3)) << std::endl; pthread_mutex_unlock(&mythreadlock); - } std::vector clusHits; @@ -781,17 +739,16 @@ namespace calc_cluster_parameter(clusHits, *my_data, iterKey->second.second); remove_hits(clusHits, rtree, adcMap); - } } void *ProcessModule(void *threadarg) { - auto my_data = static_cast(threadarg); + auto *my_data = static_cast(threadarg); ProcessModuleData(my_data); pthread_exit(nullptr); } -} //namespace +} // namespace LaserClusterizer::LaserClusterizer(const std::string &name) : SubsysReco(name) @@ -818,7 +775,7 @@ int LaserClusterizer::InitRun(PHCompositeNode *topNode) { laserClusterNodeName = "LAMINATION_CLUSTER"; } - auto laserclusters = findNode::getClass(dstNode, laserClusterNodeName); + auto *laserclusters = findNode::getClass(dstNode, laserClusterNodeName); if (!laserclusters) { PHNodeIterator dstiter(dstNode); @@ -835,7 +792,7 @@ int LaserClusterizer::InitRun(PHCompositeNode *topNode) new PHIODataNode(laserclusters, laserClusterNodeName, "PHObject"); DetNode->addNode(LaserClusterContainerNode); } - + m_geom_container = findNode::getClass(topNode, "TPCGEOMCONTAINER"); if (!m_geom_container) @@ -873,7 +830,7 @@ int LaserClusterizer::process_event(PHCompositeNode *topNode) return Fun4AllReturnCodes::ABORTRUN; } - if((eventHeader->get_RunNumber() > 66153 && !m_laserEventInfo->isGl1LaserEvent()) || (eventHeader->get_RunNumber() <= 66153 && !m_laserEventInfo->isLaserEvent())) + if ((eventHeader->get_RunNumber() > 66153 && !m_laserEventInfo->isGl1LaserEvent()) || (eventHeader->get_RunNumber() <= 66153 && !m_laserEventInfo->isLaserEvent())) { return Fun4AllReturnCodes::EVENT_OK; } @@ -897,7 +854,7 @@ int LaserClusterizer::process_event(PHCompositeNode *topNode) std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl; return Fun4AllReturnCodes::ABORTRUN; } - + // get node for clusters std::string laserClusterNodeName = "LASER_CLUSTER"; if (m_lamination) @@ -921,8 +878,8 @@ int LaserClusterizer::process_event(PHCompositeNode *topNode) return Fun4AllReturnCodes::ABORTRUN; } - TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);; - + TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId); + struct thread_pair_t { pthread_t thread{}; @@ -938,123 +895,122 @@ int LaserClusterizer::process_event(PHCompositeNode *topNode) if (pthread_mutex_init(&mythreadlock, nullptr) != 0) { - std::cout << std::endl << " mutex init failed" << std::endl; + std::cout << std::endl + << " mutex init failed" << std::endl; return 1; } - - for (unsigned int sec=0; sec<12; sec++) + + for (unsigned int sec = 0; sec < 12; sec++) { - for (int s=0; s<2; s++) + for (int s = 0; s < 2; s++) { - for (unsigned int mod=0; mod<3; mod++) + for (unsigned int mod = 0; mod < 3; mod++) { - - if(Verbosity() > 2) + if (Verbosity() > 2) { std::cout << "making thread for side: " << s << " sector: " << sec << " module: " << mod << std::endl; } - thread_pair_t &thread_pair = threads.emplace_back(); - - std::vector hitsets; - std::vector layers; - - std::vector cluster_vector; - std::vector cluster_key_vector; - - for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; - hitsetitr != hitsetrange.second; - ++hitsetitr) - { - unsigned int layer = TrkrDefs::getLayer(hitsetitr->first); - int side = TpcDefs::getSide(hitsetitr->first); - unsigned int sector = TpcDefs::getSectorId(hitsetitr->first); - if (sector != sec || side != s) - { - continue; - } - if ((mod==0 && (layer<7 || layer>22)) || (mod==1 && (layer<=22 || layer>38) ) || (mod==2 && (layer<=38 || layer>54))) - { - continue; - } - - TrkrHitSet *hitset = hitsetitr->second; - - hitsets.push_back(hitset); - layers.push_back(layer); - - } - - thread_pair.data.geom_container = m_geom_container; - thread_pair.data.tGeometry = m_tGeometry; - thread_pair.data.hitsets = hitsets; - thread_pair.data.layers = layers; - thread_pair.data.side = (bool)s; - thread_pair.data.sector = sec; + thread_pair_t &thread_pair = threads.emplace_back(); + + std::vector hitsets; + std::vector layers; + + std::vector cluster_vector; + std::vector cluster_key_vector; + + for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; + hitsetitr != hitsetrange.second; + ++hitsetitr) + { + unsigned int layer = TrkrDefs::getLayer(hitsetitr->first); + int side = TpcDefs::getSide(hitsetitr->first); + unsigned int sector = TpcDefs::getSectorId(hitsetitr->first); + if (sector != sec || side != s) + { + continue; + } + if ((mod == 0 && (layer < 7 || layer > 22)) || (mod == 1 && (layer <= 22 || layer > 38)) || (mod == 2 && (layer <= 38 || layer > 54))) + { + continue; + } + + TrkrHitSet *hitset = hitsetitr->second; + + hitsets.push_back(hitset); + layers.push_back(layer); + } + + thread_pair.data.geom_container = m_geom_container; + thread_pair.data.tGeometry = m_tGeometry; + thread_pair.data.hitsets = hitsets; + thread_pair.data.layers = layers; + thread_pair.data.side = (bool) s; + thread_pair.data.sector = sec; thread_pair.data.module = mod; - thread_pair.data.cluster_vector = cluster_vector; - thread_pair.data.cluster_key_vector = cluster_key_vector; - thread_pair.data.adc_threshold = m_adc_threshold; - thread_pair.data.peakTimeBin = m_laserEventInfo->getPeakSample(s); - thread_pair.data.layerMin = 3; - thread_pair.data.layerMax = 3; - thread_pair.data.tdriftmax = m_tdriftmax; + thread_pair.data.cluster_vector = cluster_vector; + thread_pair.data.cluster_key_vector = cluster_key_vector; + thread_pair.data.adc_threshold = m_adc_threshold; + thread_pair.data.peakTimeBin = m_laserEventInfo->getPeakSample(s); + thread_pair.data.layerMin = 3; + thread_pair.data.layerMax = 3; + thread_pair.data.tdriftmax = m_tdriftmax; thread_pair.data.eventNum = m_event; thread_pair.data.Verbosity = Verbosity(); thread_pair.data.hitHist = nullptr; thread_pair.data.doFitting = m_do_fitting; - int rc; - rc = pthread_create(&thread_pair.thread, &attr, ProcessModule, (void *) &thread_pair.data); - - if (rc) - { - std::cout << "Error:unable to create thread," << rc << std::endl; - } - - if (m_do_sequential) - { - //wait for termination of thread - int rc2 = pthread_join(thread_pair.thread, nullptr); - if (rc2) - { - std::cout << "Error:unable to join," << rc2 << std::endl; - } - - //add clusters from thread to laserClusterContainer - const auto &data(thread_pair.data); - for(int index = 0; index < (int) data.cluster_vector.size(); ++index) - { - auto cluster = data.cluster_vector[index]; - const auto ckey = data.cluster_key_vector[index]; - - m_clusterlist->addClusterSpecifyKey(ckey, cluster); - } - } + int rc; + rc = pthread_create(&thread_pair.thread, &attr, ProcessModule, (void *) &thread_pair.data); + + if (rc) + { + std::cout << "Error:unable to create thread," << rc << std::endl; + } + + if (m_do_sequential) + { + // wait for termination of thread + int rc2 = pthread_join(thread_pair.thread, nullptr); + if (rc2) + { + std::cout << "Error:unable to join," << rc2 << std::endl; + } + + // add clusters from thread to laserClusterContainer + const auto &data(thread_pair.data); + for (int index = 0; index < (int) data.cluster_vector.size(); ++index) + { + auto *cluster = data.cluster_vector[index]; + const auto ckey = data.cluster_key_vector[index]; + + m_clusterlist->addClusterSpecifyKey(ckey, cluster); + } + } } } } - + pthread_attr_destroy(&attr); if (!m_do_sequential) { - for (const auto & thread_pair : threads) + for (const auto &thread_pair : threads) { int rc2 = pthread_join(thread_pair.thread, nullptr); if (rc2) { - std::cout << "Error:unable to join," << rc2 << std::endl; + std::cout << "Error:unable to join," << rc2 << std::endl; } - - //const auto &data(thread_pair.data); - - for(int index = 0; index < (int) thread_pair.data.cluster_vector.size(); ++index) + + // const auto &data(thread_pair.data); + + for (int index = 0; index < (int) thread_pair.data.cluster_vector.size(); ++index) { - auto cluster = thread_pair.data.cluster_vector[index]; - const auto ckey = thread_pair.data.cluster_key_vector[index]; - - m_clusterlist->addClusterSpecifyKey(ckey, cluster); + auto *cluster = thread_pair.data.cluster_vector[index]; + const auto ckey = thread_pair.data.cluster_key_vector[index]; + + m_clusterlist->addClusterSpecifyKey(ckey, cluster); } } } @@ -1067,5 +1023,4 @@ int LaserClusterizer::process_event(PHCompositeNode *topNode) } return Fun4AllReturnCodes::EVENT_OK; - } diff --git a/offline/packages/tpc/LaserEventIdentifier.cc b/offline/packages/tpc/LaserEventIdentifier.cc index 8b6035db05..5968d2f09f 100644 --- a/offline/packages/tpc/LaserEventIdentifier.cc +++ b/offline/packages/tpc/LaserEventIdentifier.cc @@ -130,7 +130,7 @@ int LaserEventIdentifier::process_event(PHCompositeNode *topNode) } else if(m_runnumber > 66153) { - if ((gl1pkt->getGTMAllBusyVector() & (1<<14)) == 0) + if ((gl1pkt->getGTMAllBusyVector() & (1U<<14U)) == 0) { m_laserEventInfo->setIsGl1LaserEvent(true); m_laserEventInfo->setIsGl1LaserPileupEvent(false); diff --git a/offline/packages/tpc/Tpc3DClusterizer.cc b/offline/packages/tpc/Tpc3DClusterizer.cc index e19cbf507f..562231b71b 100644 --- a/offline/packages/tpc/Tpc3DClusterizer.cc +++ b/offline/packages/tpc/Tpc3DClusterizer.cc @@ -32,6 +32,7 @@ #include #include +#include #include #include #include @@ -72,7 +73,7 @@ int Tpc3DClusterizer::InitRun(PHCompositeNode *topNode) } // Create the Cluster node if required - auto laserclusters = findNode::getClass(dstNode, "LASER_CLUSTER"); + auto *laserclusters = findNode::getClass(dstNode, "LASER_CLUSTER"); if (!laserclusters) { PHNodeIterator dstiter(dstNode); @@ -111,13 +112,14 @@ int Tpc3DClusterizer::InitRun(PHCompositeNode *topNode) m_clusterTree->Branch("time_erase", &time_erase); m_clusterTree->Branch("time_all", &time_all); } - - if (m_output){ + + if (m_output) + { m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE"); - m_clusterNT = new TNtuple("clus3D", "clus3D","event:seed:x:y:z:r:phi:phibin:tbin:adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize"); + m_clusterNT = new TNtuple("clus3D", "clus3D", "event:seed:x:y:z:r:phi:phibin:tbin:adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize"); } - + m_geom_container = findNode::getClass(topNode, "TPCGEOMCONTAINER"); if (!m_geom_container) @@ -145,10 +147,13 @@ int Tpc3DClusterizer::InitRun(PHCompositeNode *topNode) int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) { ++m_event; - recoConsts* rc = recoConsts::instance(); - if (rc->FlagExist("RANDOMSEED")){ - m_seed = (int)rc->get_IntFlag("RANDOMSEED"); - } else { + recoConsts *rc = recoConsts::instance(); + if (rc->FlagExist("RANDOMSEED")) + { + m_seed = rc->get_IntFlag("RANDOMSEED"); + } + else + { m_seed = std::numeric_limits::quiet_NaN(); } @@ -163,11 +168,12 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) } // get node containing the digitized hits m_hits = findNode::getClass(topNode, "TRKR_HITSET"); - if (!m_hits){ + if (!m_hits) + { std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl; return Fun4AllReturnCodes::ABORTRUN; } - + // get node for clusters m_clusterlist = findNode::getClass(topNode, "LASER_CLUSTER"); if (!m_clusterlist) @@ -193,7 +199,8 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) bgi::rtree> rtree_reject; for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; - ++hitsetitr){ + ++hitsetitr) + { TrkrHitSet *hitset = hitsetitr->second; unsigned int layer = TrkrDefs::getLayer(hitsetitr->first); int side = TpcDefs::getSide(hitsetitr->first); @@ -202,49 +209,54 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) TrkrHitSet::ConstRange hitrangei = hitset->getHits(); for (TrkrHitSet::ConstIterator hitr = hitrangei.first; - hitr != hitrangei.second; - ++hitr){ + hitr != hitrangei.second; + ++hitr) + { int iphi = TpcDefs::getPad(hitr->first); int it = TpcDefs::getTBin(hitr->first); // std::cout << " iphi: " << iphi << " it: " << it << std::endl; - float_t fadc = (hitr->second->getAdc());// - m_pedestal; // proper int rounding +0.5 + float_t fadc = (hitr->second->getAdc()); // - m_pedestal; // proper int rounding +0.5 unsigned short adc = 0; - if (fadc > 0){ - adc = (unsigned short) fadc; + if (fadc > 0) + { + adc = (unsigned short) fadc; } - if (adc <= 0){ - continue; + if (adc <= 0) + { + continue; } - + std::vector testduplicate; rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001), - point(layer + 0.001, iphi + 0.001, it + 0.001))), - std::back_inserter(testduplicate)); - if (!testduplicate.empty()){ - testduplicate.clear(); - continue; + point(layer + 0.001, iphi + 0.001, it + 0.001))), + std::back_inserter(testduplicate)); + if (!testduplicate.empty()) + { + testduplicate.clear(); + continue; } TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it); - + auto spechitkey = std::make_pair(hitKey, hitsetKey); rtree_reject.insert(std::make_pair(point(1.0 * layer, 1.0 * iphi, 1.0 * it), spechitkey)); } } - + std::multimap, std::array>> adcMap; // std::cout << "n hitsets: " << std::distance(hitsetrange.first,hitsetrange.second) // << std::endl; for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; - ++hitsetitr){ + ++hitsetitr) + { TrkrHitSet *hitset = hitsetitr->second; unsigned int layer = TrkrDefs::getLayer(hitsetitr->first); int side = TpcDefs::getSide(hitsetitr->first); unsigned int sector = TpcDefs::getSectorId(hitsetitr->first); - //PHG4TpcGeom *layergeom = m_geom_container->GetLayerCellGeom(layer); - // double r = layergeom->get_radius(); - + // PHG4TpcGeom *layergeom = m_geom_container->GetLayerCellGeom(layer); + // double r = layergeom->get_radius(); + TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, side); TrkrHitSet::ConstRange hitrangei = hitset->getHits(); @@ -252,74 +264,122 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) // << std::endl; // int nhits = 0; for (TrkrHitSet::ConstIterator hitr = hitrangei.first; - hitr != hitrangei.second; - ++hitr){ + hitr != hitrangei.second; + ++hitr) + { int iphi = TpcDefs::getPad(hitr->first); int it = TpcDefs::getTBin(hitr->first); // std::cout << " iphi: " << iphi << " it: " << it << std::endl; - float_t fadc = (hitr->second->getAdc());// - m_pedestal; // proper int rounding +0.5 + float_t fadc = (hitr->second->getAdc()); // - m_pedestal; // proper int rounding +0.5 unsigned short adc = 0; // std::cout << " nhit: " << nhits++ << "adc: " << fadc << " phi: " << iphi << " it: " << it << std::endl; - if (fadc > 0){ - adc = (unsigned short) fadc; + if (fadc > 0) + { + adc = (unsigned short) fadc; } - if (adc <= 0){ - continue; + if (adc <= 0) + { + continue; } - if(layer>=7+32){ - //if(side==1)continue; - if(abs(iphi-0)<=2) continue; - if(abs(iphi-191)<=2) continue; - if(abs(iphi-206)<=1) continue; - if(abs(iphi-383)<=2) continue; - if(abs(iphi-576)<=2) continue; - if(abs(iphi-767)<=2) continue; - if(abs(iphi-960)<=2) continue; - if(abs(iphi-1522)<=2) continue; - if(abs(iphi-1344)<=2) continue; - if(abs(iphi-1536)<=2) continue; - if(abs(iphi-1728)<=2) continue; - if(abs(iphi-1920)<=2) continue; - if(abs(iphi-2111)<=2) continue; - if(abs(iphi-2303)<=2) continue; + if (layer >= 7 + 32) + { + // if(side==1)continue; + if (abs(iphi - 0) <= 2) + { + continue; + } + if (abs(iphi - 191) <= 2) + { + continue; + } + if (abs(iphi - 206) <= 1) + { + continue; + } + if (abs(iphi - 383) <= 2) + { + continue; + } + if (abs(iphi - 576) <= 2) + { + continue; + } + if (abs(iphi - 767) <= 2) + { + continue; + } + if (abs(iphi - 960) <= 2) + { + continue; + } + if (abs(iphi - 1522) <= 2) + { + continue; + } + if (abs(iphi - 1344) <= 2) + { + continue; + } + if (abs(iphi - 1536) <= 2) + { + continue; + } + if (abs(iphi - 1728) <= 2) + { + continue; + } + if (abs(iphi - 1920) <= 2) + { + continue; + } + if (abs(iphi - 2111) <= 2) + { + continue; + } + if (abs(iphi - 2303) <= 2) + { + continue; + } } /* double phi = layergeom->get_phi(iphi); double m_sampa_tbias = 39.6; double zdriftlength = (layergeom->get_zcenter(it)+ m_sampa_tbias) * m_tGeometry->get_drift_velocity(); - + float x = r * cos(phi); float y = r * sin(phi); float z = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength; if (side == 0){ - z = -z; - it = -it; + z = -z; + it = -it; } */ std::array coords = {(int) layer, iphi, it}; - + std::vector testduplicate; rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001), - point(layer + 0.001, iphi + 0.001, it + 0.001))), - std::back_inserter(testduplicate)); - if (!testduplicate.empty()){ - testduplicate.clear(); - continue; + point(layer + 0.001, iphi + 0.001, it + 0.001))), + std::back_inserter(testduplicate)); + if (!testduplicate.empty()) + { + testduplicate.clear(); + continue; } - - //test for isolated hit + + // test for isolated hit std::vector testisolated; rtree_reject.query(bgi::intersects(box(point(layer - 1.001, iphi - 1.001, it - 1.001), - point(layer + 1.001, iphi + 1.001, it + 1.001))), - std::back_inserter(testisolated)); - if(testisolated.size()==1){ - //testisolated.clear(); - continue; + point(layer + 1.001, iphi + 1.001, it + 1.001))), + std::back_inserter(testisolated)); + if (testisolated.size() == 1) + { + // testisolated.clear(); + continue; } - + TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it); - + auto spechitkey = std::make_pair(hitKey, hitsetKey); auto keyCoords = std::make_pair(spechitkey, coords); adcMap.insert(std::make_pair(adc, keyCoords)); @@ -327,38 +387,43 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) rtree.insert(std::make_pair(point(1.0 * layer, 1.0 * iphi, 1.0 * it), spechitkey)); } } - - if (Verbosity() > 1){ + + if (Verbosity() > 1) + { std::cout << "finished looping over hits" << std::endl; std::cout << "map size: " << adcMap.size() << std::endl; std::cout << "rtree size: " << rtree.size() << std::endl; } - + // done filling rTree - + t_all->restart(); - - while (adcMap.size() > 0){ + + while (!adcMap.empty()) + { auto iterKey = adcMap.rbegin(); - if (iterKey == adcMap.rend()){ + if (iterKey == adcMap.rend()) + { break; } - + auto coords = iterKey->second.second; int layer = coords[0]; int iphi = coords[1]; int it = coords[2]; - + int layerMax = layer + 1; - if (layer == 22 || layer == 38 || layer == 54){ + if (layer == 22 || layer == 38 || layer == 54) + { layerMax = layer; } int layerMin = layer - 1; - if (layer == 7 || layer == 23 || layer == 39){ + if (layer == 7 || layer == 23 || layer == 39) + { layerMin = layer; } - + std::vector clusHits; t_search->restart(); @@ -376,21 +441,24 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) clusHits.clear(); } - if (m_debug){ + if (m_debug) + { m_nClus = (int) m_eventClusters.size(); } t_all->stop(); - if (m_debug){ + if (m_debug) + { time_search = t_search->get_accumulated_time() / 1000.; time_clus = t_clus->get_accumulated_time() / 1000.; time_erase = t_erase->get_accumulated_time() / 1000.; time_all = t_all->get_accumulated_time() / 1000.; - + m_clusterTree->Fill(); } - - if (Verbosity()){ + + if (Verbosity()) + { std::cout << "rtree search time: " << t_search->get_accumulated_time() / 1000. << " sec" << std::endl; std::cout << "clustering time: " << t_clus->get_accumulated_time() / 1000. << " sec" << std::endl; std::cout << "erasing time: " << t_erase->get_accumulated_time() / 1000. << " sec" << std::endl; @@ -401,30 +469,33 @@ int Tpc3DClusterizer::process_event(PHCompositeNode *topNode) return Fun4AllReturnCodes::EVENT_OK; } -int Tpc3DClusterizer::ResetEvent(PHCompositeNode * /*topNode*/){ +int Tpc3DClusterizer::ResetEvent(PHCompositeNode * /*topNode*/) +{ m_itHist_0->Reset(); m_itHist_1->Reset(); - + if (m_debug) - { - m_tHist_0->Reset(); - m_tHist_1->Reset(); - - m_eventClusters.clear(); - } - + { + m_tHist_0->Reset(); + m_tHist_1->Reset(); + + m_eventClusters.clear(); + } + return Fun4AllReturnCodes::EVENT_OK; } int Tpc3DClusterizer::End(PHCompositeNode * /*topNode*/) { - if (m_debug){ + if (m_debug) + { m_debugFile->cd(); m_clusterTree->Write(); m_debugFile->Close(); } - if (m_output){ + if (m_output) + { m_outputFile->cd(); m_clusterNT->Write(); m_outputFile->Close(); @@ -434,7 +505,7 @@ int Tpc3DClusterizer::End(PHCompositeNode * /*topNode*/) void Tpc3DClusterizer::calc_cluster_parameter(std::vector &clusHits, std::multimap, std::array>> &adcMap) { - //std::cout << "nu clus" << std::endl; + // std::cout << "nu clus" << std::endl; double rSum = 0.0; double phiSum = 0.0; double tSum = 0.0; @@ -449,9 +520,12 @@ void Tpc3DClusterizer::calc_cluster_parameter(std::vector &clusHi TrkrDefs::hitsetkey maxKey = 0; unsigned int nHits = clusHits.size(); - int iphimin = 6666, iphimax = -1; - int ilaymin = 6666, ilaymax = -1; - float itmin = 66666666.6, itmax = -6666666666.6; + int iphimin = 6666; + int iphimax = -1; + int ilaymin = 6666; + int ilaymax = -1; + float itmin = 66666666.6; + float itmax = -6666666666.6; auto *clus = new LaserClusterv1; for (auto &clusHit : clusHits) @@ -468,21 +542,21 @@ void Tpc3DClusterizer::calc_cluster_parameter(std::vector &clusHi double phi = layergeom->get_phi(coords[1], side); double t = layergeom->get_zcenter(fabs(coords[2])); int tbin = coords[2]; - int lay = coords[0];//TrkrDefs::getLayer(spechitkey.second); + int lay = coords[0]; // TrkrDefs::getLayer(spechitkey.second); double hitzdriftlength = t * m_tGeometry->get_drift_velocity(); double hitZ = m_tdriftmax * m_tGeometry->get_drift_velocity() - hitzdriftlength; /*std::cout << " lay: " << lay - << " phi: " << phi - << " t: " << t - << " side: " << side - << std::endl; + << " phi: " << phi + << " t: " << t + << " side: " << side + << std::endl; */ - if(phiiphimax){iphimax = phi;} - if(layilaymax){ilaymax = lay;} - if(tbinitmax){itmax = tbin;} + iphimin = std::min(phi, iphimin); + iphimax = std::max(phi, iphimax); + ilaymin = std::min(lay, ilaymin); + ilaymax = std::max(lay, ilaymax); + itmin = std::min(tbin, itmin); + itmax = std::max(tbin, itmax); for (auto &iterKey : adcMap) { @@ -522,7 +596,7 @@ void Tpc3DClusterizer::calc_cluster_parameter(std::vector &clusHi if (nHits == 0) { - std::cout << "no hits"<< std::endl; + std::cout << "no hits" << std::endl; return; } @@ -554,55 +628,57 @@ void Tpc3DClusterizer::calc_cluster_parameter(std::vector &clusHi clus->setLayer(layerSum / adcSum); clus->setIPhi(iphiSum / adcSum); clus->setIT(itSum / adcSum); - int phisize = iphimax - iphimin + 1; - int lsize = ilaymax - ilaymin + 1; - int tsize = itmax - itmin +1; + int phisize = iphimax - iphimin + 1; + int lsize = ilaymax - ilaymin + 1; + int tsize = itmax - itmin + 1; if (m_debug) { m_currentCluster = (LaserCluster *) clus->CloneMe(); m_eventClusters.push_back((LaserCluster *) m_currentCluster->CloneMe()); } // if(nHits>1&&tsize>5){ - if(nHits>=1){ + if (nHits >= 1) + { const auto ckey = TrkrDefs::genClusKey(maxKey, m_clusterlist->size()); m_clusterlist->addClusterSpecifyKey(ckey, clus); - } else { + } + else + { delete clus; } - - //event:seed:x:y:z:r:phi:phibin:tbin::adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize + // event:seed:x:y:z:r:phi:phibin:tbin::adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize //"event:seed:x:y:z:r:phi:phibin:tbin::adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize" /* std::cout << " l size: " << lsize - << " phisize : " << phisize - << " tsize: " << tsize - << " maxside: " << maxside - << std::endl; + << " phisize : " << phisize + << " tsize: " << tsize + << " maxside: " << maxside + << std::endl; */ // if (m_output){ - float fX[20] = {0}; - int n = 0; - fX[n++] = m_event; - fX[n++] = m_seed; - fX[n++] = clusX; - fX[n++] = clusY; - fX[n++] = clusZ; - fX[n++] = clusR; - fX[n++] = clusPhi; - fX[n++] = clusiPhi; - fX[n++] = clusT; - fX[n++] = adcSum; - fX[n++] = maxAdc; - fX[n++] = (layerSum/adcSum); - fX[n++] = maxsector; - fX[n++] = maxside; - fX[n++] = nHits; - fX[n++] = phisize; - fX[n++] = tsize; - fX[n++] = lsize; - m_clusterNT->Fill(fX); - // } + float fX[20] = {0}; + int n = 0; + fX[n++] = m_event; + fX[n++] = m_seed; + fX[n++] = clusX; + fX[n++] = clusY; + fX[n++] = clusZ; + fX[n++] = clusR; + fX[n++] = clusPhi; + fX[n++] = clusiPhi; + fX[n++] = clusT; + fX[n++] = adcSum; + fX[n++] = maxAdc; + fX[n++] = (layerSum / adcSum); + fX[n++] = maxsector; + fX[n++] = maxside; + fX[n++] = nHits; + fX[n++] = phisize; + fX[n++] = tsize; + fX[n++] = lsize; + m_clusterNT->Fill(fX); + // } } void Tpc3DClusterizer::remove_hits(std::vector &clusHits, bgi::rtree> &rtree, std::multimap, std::array>> &adcMap) @@ -611,10 +687,11 @@ void Tpc3DClusterizer::remove_hits(std::vector &clusHits, bgi::rt { auto spechitkey = clusHit.second; - if(rtree.size()==0){ + if (rtree.empty()) + { std::cout << "not good" << std::endl; } - //rtree.remove(clusHit); + // rtree.remove(clusHit); for (auto iterAdc = adcMap.begin(); iterAdc != adcMap.end();) { @@ -623,10 +700,8 @@ void Tpc3DClusterizer::remove_hits(std::vector &clusHits, bgi::rt iterAdc = adcMap.erase(iterAdc); break; } - else - { - ++iterAdc; - } + + ++iterAdc; } } } diff --git a/offline/packages/tpc/TpcClusterMover.cc b/offline/packages/tpc/TpcClusterMover.cc index 643fe5587f..39c5e156d4 100644 --- a/offline/packages/tpc/TpcClusterMover.cc +++ b/offline/packages/tpc/TpcClusterMover.cc @@ -17,19 +17,20 @@ namespace { - [[maybe_unused]] std::ostream& operator << (std::ostream& out, const Acts::Vector3& v ) + [[maybe_unused]] std::ostream& operator<<(std::ostream& out, const Acts::Vector3& v) { out << "(" << v.x() << ", " << v.y() << ", " << v.z() << ")"; return out; } -} +} // namespace TpcClusterMover::TpcClusterMover() + : inner_tpc_spacing((mid_tpc_min_radius - inner_tpc_min_radius) / 16.0) + , mid_tpc_spacing((outer_tpc_min_radius - mid_tpc_min_radius) / 16.0) + , outer_tpc_spacing((outer_tpc_max_radius - outer_tpc_min_radius) / 16.0) { // initialize layer radii - inner_tpc_spacing = (mid_tpc_min_radius - inner_tpc_min_radius) / 16.0; - mid_tpc_spacing = (outer_tpc_min_radius - mid_tpc_min_radius) / 16.0; - outer_tpc_spacing = (outer_tpc_max_radius - outer_tpc_min_radius) / 16.0; + for (int i = 0; i < 16; ++i) { layer_radius[i] = inner_tpc_min_radius + (double) i * inner_tpc_spacing + 0.5 * inner_tpc_spacing; @@ -44,7 +45,7 @@ TpcClusterMover::TpcClusterMover() } } -void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer *cellgeo) +void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo) { if (_verbosity > 0) { @@ -65,7 +66,6 @@ void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer *cellgeo) //____________________________________________________________________________.. std::vector> TpcClusterMover::processTrack(const std::vector>& global_in) { - // Get the global positions of the TPC clusters for this track, already corrected for distortions, and move them to the surfaces // The input object contains all clusters for the track @@ -74,7 +74,7 @@ std::vector> TpcClusterMover::proces std::vector tpc_global_vec; std::vector tpc_cluskey_vec; - for (const auto& [ckey,global]:global_in) + for (const auto& [ckey, global] : global_in) { const auto trkrid = TrkrDefs::getTrkrId(ckey); if (trkrid == TrkrDefs::tpcId) @@ -85,7 +85,7 @@ std::vector> TpcClusterMover::proces else { // si clusters stay where they are - global_moved.emplace_back(ckey,global); + global_moved.emplace_back(ckey, global); } } @@ -158,7 +158,7 @@ std::vector> TpcClusterMover::proces return global_moved; } -int TpcClusterMover::get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y) +int TpcClusterMover::get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double& x, double& y) const { // finds the intersection of the fitted circle with the cylinder having radius = target_radius const auto [xplus, yplus, xminus, yminus] = TrackFitUtils::circle_circle_intersection(target_radius, R, X0, Y0); diff --git a/offline/packages/tpc/TpcClusterMover.h b/offline/packages/tpc/TpcClusterMover.h index dc67312f7f..92a053e990 100644 --- a/offline/packages/tpc/TpcClusterMover.h +++ b/offline/packages/tpc/TpcClusterMover.h @@ -27,7 +27,7 @@ class TpcClusterMover void initialize_geometry(PHG4TpcGeomContainer *cellgeo); private: - int get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y); + int get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y) const; double _z_start = 0.0; double _y_start = 0.0; diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index bd5fc3c340..c0137613f5 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -25,7 +25,6 @@ #include #include #include -#include #include #include // for SubsysReco @@ -50,6 +49,7 @@ #include +#include #include #include // for sqrt, cos, sin #include @@ -64,7 +64,7 @@ namespace { template - inline constexpr T square(const T &x) + constexpr T square(const T &x) { return x * x; } @@ -365,20 +365,14 @@ namespace int isosum = 0; int isophimin = iphi - 1; - if (isophimin < 0) - { - isophimin = 0; - } + isophimin = std::max(isophimin, 0); int isophimax = iphi + 1; if (!(isophimax < NPhiBinsMax)) { isophimax = NPhiBinsMax - 1; } int isotmin = it - 1; - if (isotmin < 0) - { - isotmin = 0; - } + isotmin = std::max(isotmin, 0); int isotmax = it + 1; if (!(isotmax < NTBinsMax)) { @@ -537,30 +531,11 @@ namespace continue; } - if (adc > max_adc) - { - max_adc = adc; - } - - if (iphi > phibinhi) - { - phibinhi = iphi; - } - - if (iphi < phibinlo) - { - phibinlo = iphi; - } - - if (it > tbinhi) - { - tbinhi = it; - } - - if (it < tbinlo) - { - tbinlo = it; - } + max_adc = std::max(max_adc, static_cast(std::round(adc))); // preserves rounding (0.5 -> 1) + phibinhi = std::max(iphi, phibinhi); + phibinlo = std::min(iphi, phibinlo); + tbinhi = std::max(it, tbinhi); + tbinlo = std::min(it, tbinlo); // if(it==it_center){ yg_sum += adc; } // update phi sums @@ -686,7 +661,7 @@ namespace // std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl; if (sqrt(phi_err_square) > my_data.min_err_squared) { - auto clus = new TrkrClusterv5; + auto *clus = new TrkrClusterv5; // auto clus = std::make_unique(); clus_base = clus; clus->setAdc(adc_sum); @@ -738,19 +713,19 @@ namespace if (my_data.fillClusHitsVerbose && b_made_cluster) { // push the data back to - my_data.phivec_ClusHitsVerbose.push_back(std::vector>{}); - my_data.zvec_ClusHitsVerbose.push_back(std::vector>{}); + my_data.phivec_ClusHitsVerbose.emplace_back(); + my_data.zvec_ClusHitsVerbose.emplace_back(); auto &vphi = my_data.phivec_ClusHitsVerbose.back(); auto &vz = my_data.zvec_ClusHitsVerbose.back(); for (auto &entry : m_phi) { - vphi.push_back({entry.first, entry.second}); + vphi.emplace_back(entry.first, entry.second); } for (auto &entry : m_z) { - vz.push_back({entry.first, entry.second}); + vz.emplace_back(entry.first, entry.second); } } @@ -875,7 +850,7 @@ namespace } if (adc > my_data->edge_threshold) { - adcval[phibin][tbin] = (unsigned short) adc; + adcval[phibin][tbin] = adc; } } } @@ -967,7 +942,7 @@ namespace } */ // std::cout << "done filling " << std::endl; - while (all_hit_map.size() > 0) + while (!all_hit_map.empty()) { // std::cout << "all hit map size: " << all_hit_map.size() << std::endl; auto iter = all_hit_map.rbegin(); @@ -1013,22 +988,10 @@ namespace { continue; } - if (wiphi > wphibinhi) - { - wphibinhi = wiphi; - } - if (wiphi < wphibinlo) - { - wphibinlo = wiphi; - } - if (wit > wtbinhi) - { - wtbinhi = wit; - } - if (wit < wtbinlo) - { - wtbinlo = wit; - } + wphibinhi = std::max(wiphi, wphibinhi); + wphibinlo = std::min(wiphi, wphibinlo); + wtbinhi = std::max(wit, wtbinhi); + wtbinlo = std::min(wit, wtbinlo); } char wtsize = wtbinhi - wtbinlo + 1; char wphisize = wphibinhi - wphibinlo + 1; @@ -1077,7 +1040,7 @@ namespace } void *ProcessSector(void *threadarg) { - auto my_data = static_cast(threadarg); + auto *my_data = static_cast(threadarg); ProcessSectorData(my_data); pthread_exit(nullptr); } @@ -1133,7 +1096,7 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) } // Create the Cluster node if required - auto trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); + auto *trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); if (!trkrclusters) { PHNodeIterator dstiter(dstNode); @@ -1151,7 +1114,7 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) DetNode->addNode(TrkrClusterContainerNode); } - auto clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); + auto *clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); if (!clusterhitassoc) { PHNodeIterator dstiter(dstNode); @@ -1168,7 +1131,7 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) DetNode->addNode(newNode); } - auto training_container = findNode::getClass(dstNode, "TRAINING_HITSET"); + auto *training_container = findNode::getClass(dstNode, "TRAINING_HITSET"); if (!training_container) { PHNodeIterator dstiter(dstNode); @@ -1217,18 +1180,18 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) if (!mClusHitsVerbose) { PHNodeIterator dstiter(dstNode); - auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); dstNode->addNode(DetNode); } mClusHitsVerbose = new ClusHitsVerbosev1(); - auto newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject"); + auto *newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject"); DetNode->addNode(newNode); } } - auto geom = + auto *geom = findNode::getClass(topNode, "TPCGEOMCONTAINER"); if (!geom) { @@ -1239,13 +1202,13 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep(); std::cout << "FirstLayerCellGeomv1 streamer: " << std::endl; - auto *g1 = (PHG4TpcGeomv1*) geom->GetFirstLayerCellGeom(); // cast because << not in the base class + auto *g1 = static_cast (geom->GetFirstLayerCellGeom()); // cast because << not in the base class std::cout << *g1 << std::endl; std::cout << "LayerCellGeomv1 streamer for layer 24: " << std::endl; - auto *g2 = (PHG4TpcGeomv1*) geom->GetLayerCellGeom(24); // cast because << not in the base class + auto *g2 = static_cast (geom->GetLayerCellGeom(24)); // cast because << not in the base class std::cout << *g2 << std::endl; std::cout << "LayerCellGeomv1 streamer for layer 40: " << std::endl; - auto *g3 = (PHG4TpcGeomv1*) geom->GetLayerCellGeom(40); // cast because << not in the base class + auto *g3 = static_cast (geom->GetLayerCellGeom(40)); // cast because << not in the base class std::cout << *g3 << std::endl; return Fun4AllReturnCodes::EVENT_OK; @@ -1489,18 +1452,18 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) const auto ckey = TrkrDefs::genClusKey(hitsetkey, index); // get cluster - auto cluster = data.cluster_vector[index]; + auto *cluster = data.cluster_vector[index]; // insert in map m_clusterlist->addClusterSpecifyKey(ckey, cluster); if (mClusHitsVerbose) { - for (auto &hit : data.phivec_ClusHitsVerbose[index]) + for (const auto &hit : data.phivec_ClusHitsVerbose[index]) { mClusHitsVerbose->addPhiHit(hit.first, hit.second); } - for (auto &hit : data.zvec_ClusHitsVerbose[index]) + for (const auto &hit : data.zvec_ClusHitsVerbose[index]) { mClusHitsVerbose->addZHit(hit.first, hit.second); } @@ -1624,7 +1587,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) const auto ckey = TrkrDefs::genClusKey(hitsetkey, index); // get cluster - auto cluster = data.cluster_vector[index]; + auto *cluster = data.cluster_vector[index]; // insert in map m_clusterlist->addClusterSpecifyKey(ckey, cluster); @@ -1668,7 +1631,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) const auto ckey = TrkrDefs::genClusKey(hitsetkey, index); // get cluster - auto cluster = data.cluster_vector[index]; + auto *cluster = data.cluster_vector[index]; // insert in map // std::cout << "X: " << cluster->getLocalX() << "Y: " << cluster->getLocalY() << std::endl; @@ -1676,11 +1639,11 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) if (mClusHitsVerbose) { - for (auto &hit : data.phivec_ClusHitsVerbose[index]) + for (const auto &hit : data.phivec_ClusHitsVerbose[index]) { mClusHitsVerbose->addPhiHit(hit.first, (float) hit.second); } - for (auto &hit : data.zvec_ClusHitsVerbose[index]) + for (const auto &hit : data.zvec_ClusHitsVerbose[index]) { mClusHitsVerbose->addZHit(hit.first, (float) hit.second); } @@ -1698,7 +1661,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) m_clusterhitassoc->addAssoc(ckey, hkey); } - for (auto v_hit : thread_pair.data.v_hits) + for (auto *v_hit : thread_pair.data.v_hits) { if (_store_hits) { diff --git a/offline/packages/tpc/TpcCombinedRawDataUnpacker.cc b/offline/packages/tpc/TpcCombinedRawDataUnpacker.cc index c194567d58..5f276c3c3a 100644 --- a/offline/packages/tpc/TpcCombinedRawDataUnpacker.cc +++ b/offline/packages/tpc/TpcCombinedRawDataUnpacker.cc @@ -57,7 +57,7 @@ void TpcCombinedRawDataUnpacker::ReadZeroSuppressedData() { m_do_zs_emulation = true; m_do_baseline_corr = false; - auto cdb = CDBInterface::instance(); + auto *cdb = CDBInterface::instance(); std::string dir = cdb->getUrl("TPC_ZS_THRESHOLDS"); auto cdbtree = std::make_unique(dir); @@ -75,7 +75,7 @@ void TpcCombinedRawDataUnpacker::ReadZeroSuppressedData() { name.str(""); name << "R"<GetSingleFloatValue(name.str().c_str()); + m_zs_threshold[i] = cdbtree->GetSingleFloatValue(name.str()); if(Verbosity() > 1) { std::cout << "Loading ADU threshold of " << m_zs_threshold[i] << " for region " << i << std::endl; diff --git a/offline/packages/tpc/TpcCombinedRawDataUnpackerDebug.cc b/offline/packages/tpc/TpcCombinedRawDataUnpackerDebug.cc index 0bb412197f..803a822aa9 100644 --- a/offline/packages/tpc/TpcCombinedRawDataUnpackerDebug.cc +++ b/offline/packages/tpc/TpcCombinedRawDataUnpackerDebug.cc @@ -34,6 +34,7 @@ #include #include +#include #include #include // for exit #include // for exit @@ -234,14 +235,8 @@ int TpcCombinedRawDataUnpackerDebug::process_event(PHCompositeNode* topNode) TpcRawHit* tpchit = tpccont->get_hit(i); uint64_t gtm_bco = tpchit->get_gtm_bco(); - if (gtm_bco < bco_min) - { - bco_min = gtm_bco; - } - if (gtm_bco > bco_max) - { - bco_max = gtm_bco; - } + bco_min = std::min(gtm_bco, bco_min); + bco_max = std::max(gtm_bco, bco_max); int fee = tpchit->get_fee(); int channel = tpchit->get_channel(); @@ -539,7 +534,7 @@ int TpcCombinedRawDataUnpackerDebug::process_event(PHCompositeNode* topNode) for (int binx = 1; binx < hist2d->GetNbinsX(); binx++) { - double timebin = ((TAxis*) hist2d->GetXaxis())->GetBinCenter(binx); + double timebin = ( hist2d->GetXaxis())->GetBinCenter(binx); std::string histname1d = "h" + std::to_string(hiter.first) + "_" + std::to_string((int) timebin); TH1D* hist1d = hist2d->ProjectionY(histname1d.c_str(), binx, binx); float local_ped = 0; @@ -697,11 +692,8 @@ int TpcCombinedRawDataUnpackerDebug::process_event(PHCompositeNode* topNode) if ((float(adc) - pedestal_offset - corr) > (hpedwidth2 * m_ped_sig_cut)) { float nuadc = (float(adc) - corr - pedestal_offset); - if (nuadc < 0) - { - nuadc = 0; - } - hitr->second->setAdc(float(nuadc)); + nuadc = std::max(nuadc, 0); + hitr->second->setAdc(nuadc); #ifdef DEBUG // hitr->second->setAdc(10); if (tbin == 383 && layer >= 7 + 32 && fee == 21) diff --git a/offline/packages/tpc/TpcDistortionCorrection.cc b/offline/packages/tpc/TpcDistortionCorrection.cc index b759331e77..71b9ea8b2f 100644 --- a/offline/packages/tpc/TpcDistortionCorrection.cc +++ b/offline/packages/tpc/TpcDistortionCorrection.cc @@ -15,7 +15,7 @@ namespace { template - inline constexpr T square(const T& x) + constexpr T square(const T x) { return x * x; } diff --git a/offline/packages/tpc/TpcLoadDistortionCorrection.cc b/offline/packages/tpc/TpcLoadDistortionCorrection.cc index 0ce2ffa4e1..65bf156a7a 100644 --- a/offline/packages/tpc/TpcLoadDistortionCorrection.cc +++ b/offline/packages/tpc/TpcLoadDistortionCorrection.cc @@ -58,7 +58,7 @@ int TpcLoadDistortionCorrection::InitRun(PHCompositeNode* topNode) } /// Get the RUN node and check - auto runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); + auto *runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); if (!runNode) { std::cout << "TpcLoadDistortionCorrection::InitRun - RUN Node missing, quitting" << std::endl; @@ -74,17 +74,17 @@ int TpcLoadDistortionCorrection::InitRun(PHCompositeNode* topNode) } // get distortion correction object and create if not found - auto distortion_correction_object = findNode::getClass(topNode, m_node_name[i]); + auto *distortion_correction_object = findNode::getClass(topNode, m_node_name[i]); if (!distortion_correction_object) { std::cout << "TpcLoadDistortionCorrection::InitRun - creating TpcDistortionCorrectionContainer in node " << m_node_name[i] << std::endl; distortion_correction_object = new TpcDistortionCorrectionContainer; - auto node = new PHDataNode(distortion_correction_object, m_node_name[i]); + auto *node = new PHDataNode(distortion_correction_object, m_node_name[i]); runNode->addNode(node); } std::cout << "TpcLoadDistortionCorrection::InitRun - reading corrections from " << m_correction_filename[i] << std::endl; - auto distortion_tfile = TFile::Open(m_correction_filename[i].c_str()); + auto *distortion_tfile = TFile::Open(m_correction_filename[i].c_str()); if (!distortion_tfile) { std::cout << "TpcLoadDistortionCorrection::InitRun - cannot open " << m_correction_filename[i] << std::endl; diff --git a/offline/packages/tpc/TpcRawDataTree.cc b/offline/packages/tpc/TpcRawDataTree.cc index 6199fb7425..1141b5de55 100644 --- a/offline/packages/tpc/TpcRawDataTree.cc +++ b/offline/packages/tpc/TpcRawDataTree.cc @@ -61,7 +61,7 @@ int TpcRawDataTree::InitRun(PHCompositeNode * /*unused*/) m_SampleTree->Branch("nWaveormInFrame", &m_nWaveormInFrame, "nWaveormInFrame/I"); m_SampleTree->Branch("maxFEECount", &m_maxFEECount, "maxFEECount/I"); m_SampleTree->Branch("nSamples", &m_nSamples, "nSamples/I"); - m_SampleTree->Branch("adcSamples", &m_adcSamples[0], "adcSamples[nSamples]/s"); + m_SampleTree->Branch("adcSamples", m_adcSamples.data(), "adcSamples[nSamples]/s"); m_SampleTree->Branch("fee", &m_fee, "fee/I"); m_SampleTree->Branch("sampaAddress", &m_sampaAddress, "sampaAddress/I"); m_SampleTree->Branch("sampaChannel", &m_sampaChannel, "sampaChannel/I"); diff --git a/offline/packages/tpc/TpcRawWriter.cc b/offline/packages/tpc/TpcRawWriter.cc index 5d7711e2bb..41c5f26df3 100644 --- a/offline/packages/tpc/TpcRawWriter.cc +++ b/offline/packages/tpc/TpcRawWriter.cc @@ -76,7 +76,7 @@ int TpcRawWriter::InitRun(PHCompositeNode *topNode) } // Create the Cluster node if required - auto trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); + auto *trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); if (!trkrclusters) { PHNodeIterator dstiter(dstNode); @@ -94,7 +94,7 @@ int TpcRawWriter::InitRun(PHCompositeNode *topNode) DetNode->addNode(TrkrClusterContainerNode); } - auto clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); + auto *clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); if (!clusterhitassoc) { PHNodeIterator dstiter(dstNode); @@ -116,7 +116,7 @@ int TpcRawWriter::InitRun(PHCompositeNode *topNode) if (!m_rawhits) { PHNodeIterator dstiter(dstNode); - auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); @@ -124,7 +124,7 @@ int TpcRawWriter::InitRun(PHCompositeNode *topNode) } m_rawhits = new RawHitSetContainerv1; - auto newNode = new PHIODataNode(m_rawhits, "TRKR_RAWHITSET", "PHObject"); + auto *newNode = new PHIODataNode(m_rawhits, "TRKR_RAWHITSET", "PHObject"); DetNode->addNode(newNode); } diff --git a/offline/packages/tpc/TpcSimpleClusterizer.cc b/offline/packages/tpc/TpcSimpleClusterizer.cc index f439d1d682..f1fb0d5b76 100644 --- a/offline/packages/tpc/TpcSimpleClusterizer.cc +++ b/offline/packages/tpc/TpcSimpleClusterizer.cc @@ -33,6 +33,7 @@ #include +#include #include #include // for sqrt, cos, sin #include @@ -46,7 +47,7 @@ namespace { template - inline constexpr T square(const T &x) + constexpr T square(const T &x) { return x * x; } @@ -142,22 +143,10 @@ namespace int iphi = iter.second.first + my_data.phioffset; int iz = iter.second.second + my_data.zoffset; - if (iphi > phibinhi) - { - phibinhi = iphi; - } - if (iphi < phibinlo) - { - phibinlo = iphi; - } - if (iz > zbinhi) - { - zbinhi = iz; - } - if (iz < zbinlo) - { - zbinlo = iz; - } + phibinhi = std::max(iphi, phibinhi); + phibinlo = std::min(iphi, phibinlo); + zbinhi = std::max(iz, zbinhi); + zbinlo = std::min(iz, zbinlo); // update phi sums double phi_center = my_data.layergeom->get_phicenter(iphi, my_data.side); @@ -205,7 +194,7 @@ namespace clusz -= (clusz < 0) ? my_data.par0_neg : my_data.par0_pos; // create cluster and fill - auto clus = new TrkrClusterv3; + auto *clus = new TrkrClusterv3; clus->setAdc(adc_sum); /// Get the surface key to find the surface from the map @@ -280,7 +269,7 @@ namespace void *ProcessSector(void *threadarg) { - auto my_data = (struct thread_data *) threadarg; + auto *my_data = (struct thread_data *) threadarg; const auto &pedestal = my_data->pedestal; const auto &phibins = my_data->phibins; @@ -332,11 +321,11 @@ namespace all_hit_map.insert(std::make_pair(adc, thisHit)); } // adcval[phibin][zbin] = (unsigned short) adc; - adcval[phibin][zbin] = (unsigned short) adc; + adcval[phibin][zbin] = adc; } } - while (all_hit_map.size() > 0) + while (!all_hit_map.empty()) { auto iter = all_hit_map.rbegin(); if (iter == all_hit_map.rend()) @@ -413,7 +402,7 @@ int TpcSimpleClusterizer::InitRun(PHCompositeNode *topNode) } // Create the Cluster node if required - auto trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); + auto *trkrclusters = findNode::getClass(dstNode, "TRKR_CLUSTER"); if (!trkrclusters) { PHNodeIterator dstiter(dstNode); @@ -431,7 +420,7 @@ int TpcSimpleClusterizer::InitRun(PHCompositeNode *topNode) DetNode->addNode(TrkrClusterContainerNode); } - auto clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); + auto *clusterhitassoc = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); if (!clusterhitassoc) { PHNodeIterator dstiter(dstNode); @@ -614,7 +603,7 @@ int TpcSimpleClusterizer::process_event(PHCompositeNode *topNode) const auto ckey = TrkrDefs::genClusKey(hitsetkey, index); // get cluster - auto cluster = data.cluster_vector[index]; + auto *cluster = data.cluster_vector[index]; // insert in map m_clusterlist->addClusterSpecifyKey(ckey, cluster); diff --git a/offline/packages/tpc/TrainingHits.cc b/offline/packages/tpc/TrainingHits.cc index 3a70db0854..20a810b225 100644 --- a/offline/packages/tpc/TrainingHits.cc +++ b/offline/packages/tpc/TrainingHits.cc @@ -1,17 +1,17 @@ #include "TrainingHits.h" TrainingHits::TrainingHits() + : radius(0.) + , phi(0.) + , z(0.) + , phistep(0.) + , zstep(0.) + , layer(0) + , ntouch(0) + , nedge(0) + , cluskey(0) { v_adc.fill(0); - radius = 0.; - phi = 0.; - z = 0.; - phistep = 0.; - zstep = 0.; - layer = 0; - ntouch = 0; - nedge = 0; - cluskey = 0; } void TrainingHits::Reset() diff --git a/offline/packages/trackbase/AlignmentTransformation.cc b/offline/packages/trackbase/AlignmentTransformation.cc index b3eb79fcb0..2e8647ab01 100644 --- a/offline/packages/trackbase/AlignmentTransformation.cc +++ b/offline/packages/trackbase/AlignmentTransformation.cc @@ -277,6 +277,7 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) surf = surfMaps.getTpcSurface(this_hitsetkey, (unsigned int) sskey); Eigen::Vector3d localFrameTranslation(0, 0, 0); + use_module_tilt = false; if (test_layer < 4 || use_module_tilt_always) { // get the local frame translation that puts the local surface center at the tilted position after the local rotations are applied @@ -285,6 +286,9 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) double this_radius = std::sqrt(this_center[0] * this_center[0] + this_center[1] * this_center[1]); float moduleRadius = TpcModuleRadii[side][sector][this_region]; // radius of the center of the module in cm localFrameTranslation = getTpcLocalFrameTranslation(moduleRadius, this_radius, sensorAngles) * 10; // cm to mm + + // set this flag for later use + use_module_tilt = true; } Acts::Transform3 transform; @@ -417,70 +421,83 @@ Acts::Transform3 AlignmentTransformation::newMakeTransform(const Surface& surf, Acts::Transform3 transform; //! If we read the survey parameters directly, that is the full transform if (survey) - { - //! The millepede affines will just be what was read in, which was the - //! survey information. This should (in principle) be equivalent to - //! the ideal position + any misalignment - transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * mpLocalRotationAffine; - } - else - { - if (trkrid == TrkrDefs::tpcId) { - transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * actsRotationAffine * mpLocalTranslationAffine * mpLocalRotationAffine; + //! The millepede affines will just be what was read in, which was the + //! survey information. This should (in principle) be equivalent to + //! the ideal position + any misalignment + transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * mpLocalRotationAffine; } - else + else { - if(use_new_silicon_rotation_order) + // not survey. this is the normal usage + + if (trkrid == TrkrDefs::tpcId) { - transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * actsRotationAffine * mpLocalTranslationAffine * mpLocalRotationAffine; + if(use_module_tilt) + { + // use module tilt transforms with local rotation followed by local translation + transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * actsRotationAffine * mpLocalTranslationAffine * mpLocalRotationAffine; + } + else + { + // backward compatibility for old alignment params sets + transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * mpLocalRotationAffine * actsRotationAffine; + } } else { - // needed for backward compatibility to existing local rotations in MVTX - transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * mpLocalRotationAffine * actsRotationAffine; + // silicon and TPOT + if(use_new_silicon_rotation_order) + { + // use new transform order for silicon as well as TPC + transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * actsRotationAffine * mpLocalTranslationAffine * mpLocalRotationAffine; + } + else + { + // needed for backward compatibility to existing local rotation parmeter sets in silicon + transform = mpGlobalTranslationAffine * mpGlobalRotationAffine * actsTranslationAffine * mpLocalRotationAffine * actsRotationAffine; + } } } - } - + if (localVerbosity) - { - Acts::Transform3 actstransform = actsTranslationAffine * actsRotationAffine; - - std::cout << "newMakeTransform" << std::endl; - std::cout << "Input sensorAngles: " << std::endl - << sensorAngles << std::endl; - std::cout << "Input sensorAnglesGlobal: " << std::endl - << sensorAnglesGlobal << std::endl; - std::cout << "Input translation: " << std::endl - << millepedeTranslation << std::endl; - std::cout << "mpLocalRotationAffine: " << std::endl - << mpLocalRotationAffine.matrix() << std::endl; - std::cout << "mpLocalTranslationAffine: " << std::endl - << mpLocalTranslationAffine.matrix() << std::endl; - std::cout << "actsRotationAffine: " << std::endl - << actsRotationAffine.matrix() << std::endl; - std::cout << "actsTranslationAffine: " << std::endl - << actsTranslationAffine.matrix() << std::endl; - std::cout << "mpRotationGlobalAffine: " << std::endl - << mpGlobalRotationAffine.matrix() << std::endl; - std::cout << "mpTranslationGlobalAffine: " << std::endl - << mpGlobalTranslationAffine.matrix() << std::endl; - std::cout << "Overall transform: " << std::endl - << transform.matrix() << std::endl; - std::cout << "overall * idealinv " << std::endl - << (transform * actstransform.inverse()).matrix() << std::endl; - std::cout << "overall - ideal " << std::endl; - for (int test = 0; test < transform.matrix().rows(); test++) { - for (int test2 = 0; test2 < transform.matrix().cols(); test2++) - { - std::cout << transform(test, test2) - actstransform(test, test2) << ", "; - } - std::cout << std::endl; + Acts::Transform3 actstransform = actsTranslationAffine * actsRotationAffine; + + std::cout << "newMakeTransform" << std::endl; + std::cout << "Input sensorAngles: " << std::endl + << sensorAngles << std::endl; + std::cout << "Input sensorAnglesGlobal: " << std::endl + << sensorAnglesGlobal << std::endl; + std::cout << "Input translation: " << std::endl + << millepedeTranslation << std::endl; + std::cout << "mpLocalRotationAffine: " << std::endl + << mpLocalRotationAffine.matrix() << std::endl; + std::cout << "mpLocalTranslationAffine: " << std::endl + << mpLocalTranslationAffine.matrix() << std::endl; + std::cout << "actsRotationAffine: " << std::endl + << actsRotationAffine.matrix() << std::endl; + std::cout << "actsTranslationAffine: " << std::endl + << actsTranslationAffine.matrix() << std::endl; + std::cout << "mpRotationGlobalAffine: " << std::endl + << mpGlobalRotationAffine.matrix() << std::endl; + std::cout << "mpTranslationGlobalAffine: " << std::endl + << mpGlobalTranslationAffine.matrix() << std::endl; + std::cout << "Overall transform: " << std::endl + << transform.matrix() << std::endl; + std::cout << "overall * idealinv " << std::endl + << (transform * actstransform.inverse()).matrix() << std::endl; + std::cout << "overall - ideal " << std::endl; + for (int test = 0; test < transform.matrix().rows(); test++) + { + for (int test2 = 0; test2 < transform.matrix().cols(); test2++) + { + std::cout << transform(test, test2) - actstransform(test, test2) << ", "; + } + std::cout << std::endl; + } } - } - + return transform; } diff --git a/offline/packages/trackbase/AlignmentTransformation.h b/offline/packages/trackbase/AlignmentTransformation.h index 7055d6c65a..e8eb3882d1 100644 --- a/offline/packages/trackbase/AlignmentTransformation.h +++ b/offline/packages/trackbase/AlignmentTransformation.h @@ -128,7 +128,8 @@ class AlignmentTransformation bool use_new_silicon_rotation_order = false; bool use_module_tilt_always = false; - + bool use_module_tilt = false; // starts at false in all cases + bool use_intt_survey_geometry = false; Acts::Transform3 newMakeTransform(const Surface& surf, Eigen::Vector3d& millepedeTranslation, Eigen::Vector3d& sensorAngles, Eigen::Vector3d& localFrameTranslation, Eigen::Vector3d& sensorAnglesGlobal, unsigned int trkrid, bool survey); diff --git a/offline/packages/trackreco/DSTClusterPruning.cc b/offline/packages/trackreco/DSTClusterPruning.cc index 15bfe547ca..36bbe4fd0f 100644 --- a/offline/packages/trackreco/DSTClusterPruning.cc +++ b/offline/packages/trackreco/DSTClusterPruning.cc @@ -177,6 +177,39 @@ void DSTClusterPruning::prune_clusters() } return; } + if(m_pruneAllSeeds) + { + for(const auto& container : {m_tpc_track_seed_container, m_silicon_track_seed_container}) + { + for (const auto& trackseed : *container) + { + if (!trackseed) + { + std::cout << "No TrackSeed" << std::endl; + continue; + } + + for (auto key_iter = trackseed->begin_cluster_keys(); key_iter != trackseed->end_cluster_keys(); ++key_iter) + { + const auto& cluster_key = *key_iter; + auto *cluster = m_cluster_map->findCluster(cluster_key); + if (!cluster) + { + std::cout << "DSTClusterPruning::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl; + continue; + } + if (!m_reduced_cluster_map->findCluster(cluster_key)) + { + m_cluster = new TrkrClusterv5(); + m_cluster->CopyFrom(cluster); + m_reduced_cluster_map->addClusterSpecifyKey(cluster_key, m_cluster); + } + } + } + } + return; + } + for (const auto& trackseed : *m_track_seed_container) { if (!trackseed) diff --git a/offline/packages/trackreco/DSTClusterPruning.h b/offline/packages/trackreco/DSTClusterPruning.h index c9ec35ea53..5fa2ff9b4f 100644 --- a/offline/packages/trackreco/DSTClusterPruning.h +++ b/offline/packages/trackreco/DSTClusterPruning.h @@ -52,6 +52,12 @@ class DSTClusterPruning : public SubsysReco //! end of processing //int End(PHCompositeNode*) override; + //! dump all clusters on all seeds out + void pruneAllSeeds() + { + m_pruneAllSeeds = true; + } + private: //! load nodes int load_nodes(PHCompositeNode*); @@ -68,6 +74,9 @@ class DSTClusterPruning : public SubsysReco TrackSeedContainer* m_tpc_track_seed_container = nullptr; TrackSeedContainer* m_silicon_track_seed_container = nullptr; +//! set to true if you want to dump out all clusters on all silicon +//! and all tpc seeds individually + bool m_pruneAllSeeds = false; //@} // debugging helpers diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index a26ee1fa0a..2ca1426301 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -5,7 +5,6 @@ * \author Tony Frawley */ - #include "PHActsTrkFitter.h" #include "ActsPropagator.h" @@ -25,7 +24,7 @@ #include #include #include -//#include +// #include #include #include #include @@ -134,8 +133,8 @@ int PHActsTrkFitter::InitRun(PHCompositeNode* topNode) MaterialSurfaceSelector selector; if (m_fitSiliconMMs || m_directNavigation) { - m_tGeometry->geometry().tGeometry->visitSurfaces(selector,false); - //std::cout<<"selector.surfaces.size() "<geometry().tGeometry->visitSurfaces(selector, false); + // std::cout<<"selector.surfaces.size() "<(m_evalname); m_evaluator->Init(topNode); - if(m_actsEvaluator && !m_simActsEvaluator) + if (m_actsEvaluator && !m_simActsEvaluator) { m_evaluator->isData(); } @@ -182,10 +181,10 @@ int PHActsTrkFitter::InitRun(PHCompositeNode* topNode) _tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); if (!_tpccellgeo) - { - std::cout << PHWHERE << " unable to find DST node TPCGEOMCONTAINER" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; - } + { + std::cout << PHWHERE << " unable to find DST node TPCGEOMCONTAINER" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } if (Verbosity() > 1) { @@ -287,7 +286,7 @@ int PHActsTrkFitter::End(PHCompositeNode* /*topNode*/) { m_evaluator->End(); } - if(m_useOutlierFinder) + if (m_useOutlierFinder) { m_outlierFinder.Write(); } @@ -314,47 +313,46 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) // capture the input crossing value, and set crossing parameters //============================== - short silicon_crossing = SHRT_MAX; - auto siseed = m_siliconSeeds->get(siid); - if(siseed) - { - silicon_crossing = siseed->get_crossing(); - } + short silicon_crossing = SHRT_MAX; + auto *siseed = m_siliconSeeds->get(siid); + if (siseed) + { + silicon_crossing = siseed->get_crossing(); + } short crossing = silicon_crossing; short int crossing_estimate = crossing; - if(m_enable_crossing_estimate) - { - crossing_estimate = track->get_crossing_estimate(); // geometric crossing estimate from matcher - } + if (m_enable_crossing_estimate) + { + crossing_estimate = track->get_crossing_estimate(); // geometric crossing estimate from matcher + } //=============================== - // must have silicon seed with valid crossing if we are doing a SC calibration fit if (m_fitSiliconMMs) + { + if ((siid == std::numeric_limits::max()) || (silicon_crossing == SHRT_MAX)) { - if( (siid == std::numeric_limits::max()) || (silicon_crossing == SHRT_MAX)) - { - continue; - } + continue; } + } // do not skip TPC only tracks, just set crossing to the nominal zero - if(!siseed) - { - crossing = 0; - } + if (!siseed) + { + crossing = 0; + } if (Verbosity() > 1) { - if(siseed) - { - std::cout << "tpc and si id " << tpcid << ", " << siid << " silicon_crossing " << silicon_crossing - << " crossing " << crossing << " crossing estimate " << crossing_estimate << std::endl; - } + if (siseed) + { + std::cout << "tpc and si id " << tpcid << ", " << siid << " silicon_crossing " << silicon_crossing + << " crossing " << crossing << " crossing estimate " << crossing_estimate << std::endl; + } } - auto tpcseed = m_tpcSeeds->get(tpcid); + auto *tpcseed = m_tpcSeeds->get(tpcid); /// Need to also check that the tpc seed wasn't removed by the ghost finder if (!tpcseed) @@ -381,7 +379,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) if (Verbosity() > 1 && siseed) { std::cout << " m_pp_mode " << m_pp_mode << " m_enable_crossing_estimate " << m_enable_crossing_estimate - << " INTT crossing " << crossing << " crossing_estimate " << crossing_estimate << std::endl; + << " INTT crossing " << crossing << " crossing_estimate " << crossing_estimate << std::endl; } short int this_crossing = crossing; @@ -390,35 +388,35 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) std::vector chisq_ndf; std::vector svtx_vec; - if(m_pp_mode) + if (m_pp_mode) + { + if (m_enable_crossing_estimate && crossing == SHRT_MAX) + { + // this only happens if there is a silicon seed but no assigned INTT crossing, and only in pp_mode + // If there is no INTT crossing, start with the crossing_estimate value, vary up and down, fit, and choose the best chisq/ndf + use_estimate = true; + nvary = max_bunch_search; + if (Verbosity() > 1) + { + std::cout << " No INTT crossing: use crossing_estimate " << crossing_estimate << " with nvary " << nvary << std::endl; + } + } + else { - if (m_enable_crossing_estimate && crossing == SHRT_MAX) - { - // this only happens if there is a silicon seed but no assigned INTT crossing, and only in pp_mode - // If there is no INTT crossing, start with the crossing_estimate value, vary up and down, fit, and choose the best chisq/ndf - use_estimate = true; - nvary = max_bunch_search; - if (Verbosity() > 1) - { - std::cout << " No INTT crossing: use crossing_estimate " << crossing_estimate << " with nvary " << nvary << std::endl; - } - } - else - { - // use INTT crossing - crossing_estimate = crossing; - } + // use INTT crossing + crossing_estimate = crossing; } + } else + { + // non pp mode, we want only crossing zero, veto others + if (siseed && silicon_crossing != 0) { - // non pp mode, we want only crossing zero, veto others - if(siseed && silicon_crossing != 0) - { - crossing = 0; - //continue; - } - crossing_estimate = crossing; + crossing = 0; + // continue; } + crossing_estimate = crossing; + } // Fit this track assuming either: // crossing = INTT value, if it exists (uses nvary = 0) @@ -441,16 +439,16 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) makeSourceLinks.initialize(_tpccellgeo); makeSourceLinks.setVerbosity(Verbosity()); makeSourceLinks.set_pp_mode(m_pp_mode); - for(const auto& layer : m_ignoreLayer) + for (const auto& layer : m_ignoreLayer) { makeSourceLinks.ignoreLayer(layer); } // loop over modifiedTransformSet and replace transient elements modified for the previous track with the default transforms // does nothing if m_transient_id_set is empty makeSourceLinks.resetTransientTransformMap( - m_alignmentTransformationMapTransient, - m_transient_id_set, - m_tGeometry); + m_alignmentTransformationMapTransient, + m_transient_id_set, + m_tGeometry); if (m_use_clustermover) { @@ -459,37 +457,56 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) { // silicon source links sourceLinks = makeSourceLinks.getSourceLinksClusterMover( - siseed, + siseed, + measurements, + m_clusterContainer, + m_tGeometry, + m_globalPositionWrapper, + this_crossing); + } + + // tpc source links + const auto tpcSourceLinks = makeSourceLinks.getSourceLinksClusterMover( + tpcseed, measurements, m_clusterContainer, m_tGeometry, m_globalPositionWrapper, this_crossing); - } - - // tpc source links - const auto tpcSourceLinks = makeSourceLinks.getSourceLinksClusterMover( - tpcseed, - measurements, - m_clusterContainer, - m_tGeometry, - m_globalPositionWrapper, - this_crossing); // add tpc sourcelinks to silicon source links sourceLinks.insert(sourceLinks.end(), tpcSourceLinks.begin(), tpcSourceLinks.end()); - - } else { - + } + else + { // make source links using transient transforms for distortion corrections - if(Verbosity() > 1) - { std::cout << "Calling getSourceLinks for si seed, siid " << siid << " and tpcid " << tpcid << std::endl; } + if (Verbosity() > 1) + { + std::cout << "Calling getSourceLinks for si seed, siid " << siid << " and tpcid " << tpcid << std::endl; + } if (siseed && !m_ignoreSilicon) { // silicon source links sourceLinks = makeSourceLinks.getSourceLinks( - siseed, + siseed, + measurements, + m_clusterContainer, + m_tGeometry, + m_globalPositionWrapper, + m_alignmentTransformationMapTransient, + m_transient_id_set, + this_crossing); + } + + if (Verbosity() > 1) + { + std::cout << "Calling getSourceLinks for tpc seed, siid " << siid << " and tpcid " << tpcid << std::endl; + } + + // tpc source links + const auto tpcSourceLinks = makeSourceLinks.getSourceLinks( + tpcseed, measurements, m_clusterContainer, m_tGeometry, @@ -497,21 +514,6 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) m_alignmentTransformationMapTransient, m_transient_id_set, this_crossing); - } - - if(Verbosity() > 1) - { std::cout << "Calling getSourceLinks for tpc seed, siid " << siid << " and tpcid " << tpcid << std::endl; } - - // tpc source links - const auto tpcSourceLinks = makeSourceLinks.getSourceLinks( - tpcseed, - measurements, - m_clusterContainer, - m_tGeometry, - m_globalPositionWrapper, - m_alignmentTransformationMapTransient, - m_transient_id_set, - this_crossing); // add tpc sourcelinks to silicon source links sourceLinks.insert(sourceLinks.end(), tpcSourceLinks.begin(), tpcSourceLinks.end()); @@ -524,15 +526,15 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) Acts::Vector3 position(0, 0, 0); if (siseed) { - position = TrackSeedHelper::get_xyz(siseed)*Acts::UnitConstants::cm; + position = TrackSeedHelper::get_xyz(siseed) * Acts::UnitConstants::cm; } - if(!siseed || !is_valid(position) || m_ignoreSilicon) + if (!siseed || !is_valid(position) || m_ignoreSilicon) { - position = TrackSeedHelper::get_xyz(tpcseed)*Acts::UnitConstants::cm; + position = TrackSeedHelper::get_xyz(tpcseed) * Acts::UnitConstants::cm; } if (!is_valid(position)) { - if(Verbosity() > 4) + if (Verbosity() > 4) { std::cout << "Invalid position of " << position.transpose() << std::endl; } @@ -559,26 +561,26 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) for (const auto& surface_apr : m_materialSurfaces) { - if(m_forceSiOnlyFit) + if (m_forceSiOnlyFit) { - if(surface_apr->geometryId().volume() >12) + if (surface_apr->geometryId().volume() > 12) { continue; } } bool pop_flag = false; - if(surface_apr->geometryId().approach() == 1) + if (surface_apr->geometryId().approach() == 1) { surfaces.push_back(surface_apr); } else { pop_flag = true; - for (const auto& surface_sns: surfaces_tmp) + for (const auto& surface_sns : surfaces_tmp) { if (surface_apr->geometryId().volume() == surface_sns->geometryId().volume()) { - if ( surface_apr->geometryId().layer()==surface_sns->geometryId().layer()) + if (surface_apr->geometryId().layer() == surface_sns->geometryId().layer()) { pop_flag = false; surfaces.push_back(surface_sns); @@ -594,9 +596,9 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) surfaces.pop_back(); pop_flag = false; } - if (surface_apr->geometryId().volume() == 12&& surface_apr->geometryId().layer()==8) + if (surface_apr->geometryId().volume() == 12 && surface_apr->geometryId().layer() == 8) { - for (const auto& surface_sns: surfaces_tmp) + for (const auto& surface_sns : surfaces_tmp) { if (14 == surface_sns->geometryId().volume()) { @@ -619,13 +621,13 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) { // make sure micromegas are in the tracks, if required if (m_useMicromegas && - std::none_of(surfaces.begin(), surfaces.end(), [this](const auto& surface) - { return m_tGeometry->maps().isMicromegasSurface(surface); })) - { - continue; + std::none_of(surfaces.begin(), surfaces.end(), [this](const auto& surface) + { return m_tGeometry->maps().isMicromegasSurface(surface); })) + { + continue; + } } } - } float px = std::numeric_limits::quiet_NaN(); float py = std::numeric_limits::quiet_NaN(); @@ -635,7 +637,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) float seedphi = 0; float seedtheta = 0; float seedeta = 0; - if(siseed) + if (siseed) { seedphi = siseed->get_phi(); seedtheta = siseed->get_theta(); @@ -659,7 +661,9 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) px = pt * std::cos(phi); py = pt * std::sin(phi); pz = pt * std::cosh(eta) * std::cos(theta); - } else { + } + else + { px = seedpt * std::cos(seedphi); py = seedpt * std::sin(seedphi); pz = seedpt * std::cosh(seedeta) * std::cos(seedtheta); @@ -668,14 +672,14 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) Acts::Vector3 momentum(px, py, pz); if (!is_valid(momentum)) { - if(Verbosity() > 4) + if (Verbosity() > 4) { std::cout << "Invalid momentum of " << momentum.transpose() << std::endl; } continue; } - auto pSurface = Acts::Surface::makeShared( position); + auto pSurface = Acts::Surface::makeShared(position); Acts::Vector4 actsFourPos(position(0), position(1), position(2), 10 * Acts::UnitConstants::ns); Acts::BoundSquareMatrix cov = setDefaultCovariance(); @@ -723,8 +727,10 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) auto trackStateContainer = std::make_shared(); ActsTrackFittingAlgorithm::TrackContainer tracks(trackContainer, trackStateContainer); - if(Verbosity() > 1) - { std::cout << "Calling fitTrack for track with siid " << siid << " tpcid " << tpcid << " crossing " << crossing << std::endl; } + if (Verbosity() > 1) + { + std::cout << "Calling fitTrack for track with siid " << siid << " tpcid " << tpcid << " crossing " << crossing << std::endl; + } auto result = fitTrack(sourceLinks, seed, kfOptions, surfaces, calibrator, tracks); fitTimer.stop(); @@ -761,7 +767,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) if (ivary != nvary) { - if(Verbosity() > 3) + if (Verbosity() > 3) { std::cout << "Skipping track fit for trial variation" << std::endl; } @@ -806,7 +812,6 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) if (getTrackFitResult(result, track, &newTrack, tracks, measurements)) { - // insert in dedicated map m_directedTrackMap->insertWithKey(&newTrack, trid); } @@ -822,11 +827,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) m_trackMap->insertWithKey(&newTrack, trid); } } // end insert track for normal fit - } // end case where INTT crossing is known - - - - + } // end case where INTT crossing is known } else if (!m_fitSiliconMMs) { @@ -840,7 +841,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) << std::endl; } } // end fit failed case - } // end ivary loop + } // end ivary loop trackTimer.stop(); auto trackTime = trackTimer.get_accumulated_time(); @@ -855,10 +856,10 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) } bool PHActsTrkFitter::getTrackFitResult( - const FitResult& fitOutput, - TrackSeed* seed, SvtxTrack* track, - const ActsTrackFittingAlgorithm::TrackContainer& tracks, - const ActsTrackFittingAlgorithm::MeasurementContainer& measurements) + const FitResult& fitOutput, + TrackSeed* seed, SvtxTrack* track, + const ActsTrackFittingAlgorithm::TrackContainer& tracks, + const ActsTrackFittingAlgorithm::MeasurementContainer& measurements) { /// Make a trajectory state for storage, which conforms to Acts track fit /// analysis tool @@ -872,12 +873,12 @@ bool PHActsTrkFitter::getTrackFitResult( // retrieve track parameters from fit result Acts::BoundTrackParameters parameters = ActsExamples::TrackParameters(outtrack.referenceSurface().getSharedPtr(), - outtrack.parameters(), outtrack.covariance(), outtrack.particleHypothesis()); + outtrack.parameters(), outtrack.covariance(), outtrack.particleHypothesis()); indexedParams.emplace( - outtrack.tipIndex(), - ActsExamples::TrackParameters{outtrack.referenceSurface().getSharedPtr(), - outtrack.parameters(), outtrack.covariance(), outtrack.particleHypothesis()}); + outtrack.tipIndex(), + ActsExamples::TrackParameters{outtrack.referenceSurface().getSharedPtr(), + outtrack.parameters(), outtrack.covariance(), outtrack.particleHypothesis()}); if (Verbosity() > 2) { @@ -948,7 +949,9 @@ ActsTrackFittingAlgorithm::TrackFitterResult PHActsTrkFitter::fitTrack( { // use direct fit for silicon MM gits or direct navigation if (m_fitSiliconMMs || m_directNavigation) - { return (*m_fitCfg.dFit)(sourceLinks, seed, kfOptions, surfSequence, calibrator, tracks); } + { + return (*m_fitCfg.dFit)(sourceLinks, seed, kfOptions, surfSequence, calibrator, tracks); + } // use full fit in all other cases return (*m_fitCfg.fit)(sourceLinks, seed, kfOptions, calibrator, tracks); @@ -986,9 +989,9 @@ SourceLinkVec PHActsTrkFitter::getSurfaceVector(const SourceLinkVec& sourceLinks } } - if(m_forceSiOnlyFit) + if (m_forceSiOnlyFit) { - if(m_tGeometry->maps().isMicromegasSurface(surf)||m_tGeometry->maps().isTpcSurface(surf)) + if (m_tGeometry->maps().isMicromegasSurface(surf) || m_tGeometry->maps().isTpcSurface(surf)) { continue; } @@ -1059,10 +1062,10 @@ void PHActsTrkFitter::checkSurfaceVec(SurfacePtrVec& surfaces) const } void PHActsTrkFitter::updateSvtxTrack( - const std::vector& tips, - const Trajectory::IndexedParameters& paramsMap, - const ActsTrackFittingAlgorithm::TrackContainer& tracks, - SvtxTrack* track) + const std::vector& tips, + const Trajectory::IndexedParameters& paramsMap, + const ActsTrackFittingAlgorithm::TrackContainer& tracks, + SvtxTrack* track) { const auto& mj = tracks.trackStateContainer(); @@ -1133,31 +1136,37 @@ void PHActsTrkFitter::updateSvtxTrack( trackStateTimer.restart(); if (m_fillSvtxTrackStates) - { transformer.fillSvtxTrackStates(mj, trackTip, track, m_transient_geocontext); } + { + transformer.fillSvtxTrackStates(mj, trackTip, track, m_transient_geocontext); + } // in using silicon mm fit also extrapolate track parameters to all TPC surfaces with clusters // get all tpc clusters auto* seed = track->get_tpc_seed(); - if( m_fitSiliconMMs && seed ) + if (m_fitSiliconMMs && seed) { - // acts propagator ActsPropagator propagator(m_tGeometry); // loop over cluster keys associated to TPC seed - for( auto key_iter = seed->begin_cluster_keys(); key_iter != seed->end_cluster_keys(); ++key_iter ) + for (auto key_iter = seed->begin_cluster_keys(); key_iter != seed->end_cluster_keys(); ++key_iter) { const auto& cluskey = *key_iter; // make sure cluster is from TPC const auto detId = TrkrDefs::getTrkrId(cluskey); if (detId != TrkrDefs::tpcId) - { continue; } + { + continue; + } // get layer, propagate const auto layer = TrkrDefs::getLayer(cluskey); auto result = propagator.propagateTrack(params, layer); - if( !result.ok() ) { continue; } + if (!result.ok()) + { + continue; + } // get path length and extrapolated parameters auto& [pathLength, trackStateParams] = result.value(); diff --git a/offline/packages/trackreco/PHActsTrkFitter.h b/offline/packages/trackreco/PHActsTrkFitter.h index c6e0afec35..cca5095c79 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.h +++ b/offline/packages/trackreco/PHActsTrkFitter.h @@ -21,8 +21,8 @@ #include #include #include -#include #include +#include #include @@ -130,18 +130,19 @@ class PHActsTrkFitter : public SubsysReco void set_track_map_name(const std::string& map_name) { _track_map_name = map_name; } void set_svtx_seed_map_name(const std::string& map_name) { _svtx_seed_map_name = map_name; } - void set_svtx_alignment_state_map_name(const std::string& map_name) { - _svtx_alignment_state_map_name = map_name; - m_alignStates.alignmentStateMap(map_name); + void set_svtx_alignment_state_map_name(const std::string& map_name) + { + _svtx_alignment_state_map_name = map_name; + m_alignStates.alignmentStateMap(map_name); } /// Set flag for pp running void set_pp_mode(bool ispp) { m_pp_mode = ispp; } - void set_enable_geometric_crossing_estimate(bool flag) { m_enable_crossing_estimate = flag ; } + void set_enable_geometric_crossing_estimate(bool flag) { m_enable_crossing_estimate = flag; } void set_use_clustermover(bool use) { m_use_clustermover = use; } void ignoreLayer(int layer) { m_ignoreLayer.insert(layer); } - void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; } + void setTrkrClusterContainerName(const std::string& name) { m_clusterContainerName = name; } void setDirectNavigation(bool flag) { m_directNavigation = flag; } private: @@ -155,10 +156,10 @@ class PHActsTrkFitter : public SubsysReco /// Convert the acts track fit result to an svtx track void updateSvtxTrack( - const std::vector& tips, - const Trajectory::IndexedParameters& paramsMap, - const ActsTrackFittingAlgorithm::TrackContainer& tracks, - SvtxTrack* track); + const std::vector& tips, + const Trajectory::IndexedParameters& paramsMap, + const ActsTrackFittingAlgorithm::TrackContainer& tracks, + SvtxTrack* track); /// Helper function to call either the regular navigation or direct /// navigation, depending on m_fitSiliconMMs @@ -240,7 +241,7 @@ class PHActsTrkFitter : public SubsysReco // max variation of bunch crossing away from crossing_estimate short int max_bunch_search = 2; - //name of TRKR_CLUSTER container + // name of TRKR_CLUSTER container std::string m_clusterContainerName = "TRKR_CLUSTER"; //!@name evaluator @@ -253,7 +254,7 @@ class PHActsTrkFitter : public SubsysReco //@} //! tracks -// SvtxTrackMap* m_seedTracks = nullptr; + // SvtxTrackMap* m_seedTracks = nullptr; //! tpc global position wrapper TpcGlobalPositionWrapper m_globalPositionWrapper; @@ -268,7 +269,7 @@ class PHActsTrkFitter : public SubsysReco int _n_iteration = 0; std::string _track_map_name = "SvtxTrackMap"; std::string _svtx_seed_map_name = "SvtxTrackSeedContainer"; - std::string _svtx_alignment_state_map_name = "SvtxAlignmentStateMap"; + std::string _svtx_alignment_state_map_name = "SvtxAlignmentStateMap"; /// Default particle assumption to pion unsigned int m_pHypothesis = 211; @@ -292,14 +293,18 @@ class PHActsTrkFitter : public SubsysReco std::vector m_materialSurfaces = {}; - struct MaterialSurfaceSelector { + struct MaterialSurfaceSelector + { std::vector surfaces = {}; /// @param surface is the test surface - void operator()(const Acts::Surface* surface) { - if (surface->surfaceMaterial() != nullptr) { + void operator()(const Acts::Surface* surface) + { + if (surface->surfaceMaterial() != nullptr) + { if (std::find(surfaces.begin(), surfaces.end(), surface) == - surfaces.end()) { + surfaces.end()) + { surfaces.push_back(surface); } } diff --git a/offline/packages/trackreco/PHSiliconTpcTrackMatching.cc b/offline/packages/trackreco/PHSiliconTpcTrackMatching.cc index 2e66a6df28..6aaae4bf93 100644 --- a/offline/packages/trackreco/PHSiliconTpcTrackMatching.cc +++ b/offline/packages/trackreco/PHSiliconTpcTrackMatching.cc @@ -445,10 +445,10 @@ int PHSiliconTpcTrackMatching::GetNodes(PHCompositeNode *topNode) svtxNode->addNode(node); } - _cluster_map = findNode::getClass(topNode, "TRKR_CLUSTER"); + _cluster_map = findNode::getClass(topNode, _cluster_map_name); if (!_cluster_map) { - std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl; + std::cout << PHWHERE << " ERROR: Can't find node " <<_cluster_map_name << std::endl; return Fun4AllReturnCodes::ABORTEVENT; } diff --git a/offline/packages/trackreco/PHSiliconTpcTrackMatching.h b/offline/packages/trackreco/PHSiliconTpcTrackMatching.h index acb157bd83..6de4a1dd3a 100644 --- a/offline/packages/trackreco/PHSiliconTpcTrackMatching.h +++ b/offline/packages/trackreco/PHSiliconTpcTrackMatching.h @@ -139,7 +139,10 @@ class PHSiliconTpcTrackMatching : public SubsysReco, public PHParameterInterface void set_file_name(const std::string &name) { _file_name = name; } void set_pp_mode(const bool flag) { _pp_mode = flag; } void set_use_intt_crossing(const bool flag) { _use_intt_crossing = flag; } - + void set_cluster_map_name(const std::string &name) + { + _cluster_map_name = name; + } int InitRun(PHCompositeNode *topNode) override; int process_event(PHCompositeNode *) override; @@ -210,6 +213,7 @@ class PHSiliconTpcTrackMatching : public SubsysReco, public PHParameterInterface int _n_iteration = 0; std::string _track_map_name = "TpcTrackSeedContainer"; std::string _silicon_track_map_name = "SiliconTrackSeedContainer"; + std::string _cluster_map_name = "TRKR_CLUSTER"; std::string m_fieldMap = "1.4"; std::vector getTrackletClusterList(TrackSeed* tracklet); }; diff --git a/offline/packages/trackreco/PHSimpleVertexFinder.h b/offline/packages/trackreco/PHSimpleVertexFinder.h index da0e36fabf..6915693c4b 100644 --- a/offline/packages/trackreco/PHSimpleVertexFinder.h +++ b/offline/packages/trackreco/PHSimpleVertexFinder.h @@ -54,7 +54,7 @@ class PHSimpleVertexFinder : public SubsysReco void setTrackMapName(const std::string &name) { _track_map_name = name; } void setVertexMapName(const std::string &name) { _vertex_map_name = name; } void zeroField(const bool flag) { _zero_field = flag; } - void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; } + void setTrkrClusterContainerName(const std::string &name){ m_clusterContainerName = name; } void set_pp_mode(bool mode) { _pp_mode = mode; } private: diff --git a/offline/packages/trackreco/PHTpcDeltaZCorrection.h b/offline/packages/trackreco/PHTpcDeltaZCorrection.h index 46d099f587..4d383cc854 100644 --- a/offline/packages/trackreco/PHTpcDeltaZCorrection.h +++ b/offline/packages/trackreco/PHTpcDeltaZCorrection.h @@ -32,7 +32,7 @@ class PHTpcDeltaZCorrection : public SubsysReco, public PHParameterInterface int process_event(PHCompositeNode *topNode) override; int End(PHCompositeNode *topNode) override; void SetDefaultParameters() override; - void setTrkrClusterContainerName(std::string &name) { m_clusterContainerName = name; } + void setTrkrClusterContainerName(const std::string &name) { m_clusterContainerName = name; } private: /// load nodes