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