diff --git a/otsdaq-mu2e-calorimeter/ArtModules/CMakeLists.txt b/otsdaq-mu2e-calorimeter/ArtModules/CMakeLists.txt index 1539eea..4b88cb3 100755 --- a/otsdaq-mu2e-calorimeter/ArtModules/CMakeLists.txt +++ b/otsdaq-mu2e-calorimeter/ArtModules/CMakeLists.txt @@ -1,3 +1,8 @@ +include_directories($ENV{KINKAL_INC}) +include_directories($ENV{BTRK_INC}) +link_directories($ENV{KINKAL_LIB}) +link_directories($ENV{BTRK_LIB}) + cet_build_plugin(testFragments art::module LIBRARIES REG artdaq_core_mu2e::artdaq-core-mu2e_Overlays art_root_io::TFileService_service ROOT::Hist ) @@ -9,15 +14,33 @@ cet_build_plugin(CaloDataVerifier art::module LIBRARIES REG artdaq::DAQdata ) -cet_build_plugin(CaloDataAnalyzer art::module LIBRARIES REG - artdaq_core_mu2e::artdaq-core-mu2e_Data - artdaq_core_mu2e::artdaq-core-mu2e_Data_dict - artdaq_core_mu2e::artdaq-core-mu2e_Overlays - artdaq::DAQdata - art_root_io::TFileService_service - ROOT::Hist - ROOT::Tree +cet_build_plugin(CaloDataAnalyzer art::module + REG_SOURCE CaloDataAnalyzer_module.cc + LIBRARIES REG + artdaq_core_mu2e::artdaq-core-mu2e_Data + artdaq_core_mu2e::artdaq-core-mu2e_Data_dict + artdaq_core_mu2e::artdaq-core-mu2e_Overlays + artdaq::DAQdata + art_root_io::TFileService_service + Offline::DAQ + ROOT::Hist + ROOT::Tree +) + +cet_build_plugin(CaloiercAnalyzer art::module + REG_SOURCE CaloiercAnalyzer_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::RecoDataProducts + Offline::DAQ + Offline::Mu2eKinKal + ROOT::Hist + ROOT::Tree + ROOT::Core + ROOT::RIO + ROOT::Gui ) + install_headers() install_source() diff --git a/otsdaq-mu2e-calorimeter/ArtModules/CaloDataAnalyzer_module.cc b/otsdaq-mu2e-calorimeter/ArtModules/CaloDataAnalyzer_module.cc index 0478058..2e2f807 100644 --- a/otsdaq-mu2e-calorimeter/ArtModules/CaloDataAnalyzer_module.cc +++ b/otsdaq-mu2e-calorimeter/ArtModules/CaloDataAnalyzer_module.cc @@ -16,6 +16,7 @@ #include #include "artdaq-core/Data/Fragment.hh" +#include "Offline/DAQ/inc/CaloDAQUtilities.hh" #include "artdaq-core-mu2e/Data/CalorimeterDataDecoder.hh" #include "artdaq-core-mu2e/Data/EventHeader.hh" #include "artdaq-core-mu2e/Overlays/DTCEventFragment.hh" @@ -38,34 +39,36 @@ namespace mu2e class CaloDataAnalyzer : public art::EDAnalyzer { public: - struct Config - { - fhicl::Atom unpackerModuleLabel{ - fhicl::Name("unpackerModuleLabel"), - fhicl::Comment("unpackerModuleLabel"), - ""}; - fhicl::Atom unpackerInstanceLabel{ - fhicl::Name("unpackerInstanceLabel"), - fhicl::Comment("unpackerInstanceLabel"), - ""}; - fhicl::Atom verbosity{ - fhicl::Name("verbosity"), fhicl::Comment("Verbosity [0-2]"), 0}; - fhicl::Atom data_type{ - fhicl::Name("dataType"), - fhicl::Comment("Data type (0:standard, 1:debug, 2:counters)"), - 0}; - }; + // clang-format off + struct Config { + fhicl::Atom unpackerModuleLabel {fhicl::Name("unpackerModuleLabel" ) , fhicl::Comment("unpackerModuleLabel"), ""}; + fhicl::Atom unpackerInstanceLabel {fhicl::Name("unpackerInstanceLabel" ) , fhicl::Comment("unpackerInstanceLabel"), ""}; + fhicl::Atom verbosity {fhicl::Name("verbosity" ) , fhicl::Comment("Verbosity [0-2]"), 0}; + fhicl::Atom data_type {fhicl::Name("dataType" ) , fhicl::Comment("Data type (0:standard, 1:debug, 2:counters)"), 0}; + fhicl::Atom maxEventNum {fhicl::Name("maxEventNum" ) , fhicl::Comment("maxEventNum (-1:infinite)"), -1}; + fhicl::Atom fillEmptyEvents {fhicl::Name("fillEmptyEvents" ) , fhicl::Comment("Fill tree even if zero hits"), false}; + }; + // clang-format on explicit CaloDataAnalyzer(const art::EDAnalyzer::Table& config); void analyze(art::Event const& event) override; + void endJob() override; void processCaloData(mu2e::CalorimeterDataDecoder const& caloDecoder); private: - std::string unpackerModuleLabel_; - std::string unpackerInstanceLabel_; - int verbosity_; - int data_type_; - + std::string unpackerModuleLabel_; + std::string unpackerInstanceLabel_; + int verbosity_; + int data_type_; + mu2e::CaloDAQUtilities caloDAQUtil_; + int maxEventNum_; + int fillEmptyEvents_; + + int event_failedhits; + int total_events; + int total_hits; + int total_goodhits; + int total_failedhits; size_t nCaloEvents; size_t nCaloHits; int this_eventNumber; @@ -80,31 +83,29 @@ class CaloDataAnalyzer : public art::EDAnalyzer TGraph* g_eventEWT; static const int nROCs = 6; - static const int nCHs = 20; - static const int nCHANs = 2; static const int MAXNHITS = 150; static const int MAXNSAMPLES = 6300; - int t_run; - int t_subrun; - int t_nevt; - int t_dtcID; - Long64_t t_currentDTCEventWindow; - Long64_t t_currentROCEventWindow[nROCs]; - int t_nhits; - int t_boardID[MAXNHITS]; - int t_linkID[MAXNHITS]; - int t_chanID[MAXNHITS]; - int t_errflag[MAXNHITS]; - int t_fff[MAXNHITS]; - int t_time_tot[MAXNHITS]; - int t_ewhit[MAXNHITS]; - int t_peakpos[MAXNHITS]; - int t_peakval[MAXNHITS]; - int t_nofsamples[MAXNHITS]; - int t_firstsample[MAXNHITS]; - int t_nsamples; - int t_ADC[MAXNSAMPLES]; + int t_run; + int t_subrun; + int t_nevt; + Long64_t t_currentDTCEventWindow; + int t_nhits; + int t_dtcID[MAXNHITS]; + int t_boardID[MAXNHITS]; + int t_linkID[MAXNHITS]; + int t_chanID[MAXNHITS]; + int t_errflag[MAXNHITS]; + int t_fff[MAXNHITS]; + int t_time_tot[MAXNHITS]; + Long64_t t_ewhit[MAXNHITS]; + int t_peakpos[MAXNHITS]; + int t_peakval[MAXNHITS]; + int t_nofsamples[MAXNHITS]; + int t_firstsample[MAXNHITS]; + int t_nsamples; + int t_ADC[MAXNSAMPLES]; + std::vector>* t_ADC_hit = 0; int tH_run; int tH_subrun; @@ -119,7 +120,7 @@ class CaloDataAnalyzer : public art::EDAnalyzer int tH_errflag; int tH_fff; int tH_time_tot; - int tH_ewhit; + Long64_t tH_ewhit; int tH_peakpos; int tH_peakval; int tH_nofsamples; @@ -136,6 +137,9 @@ mu2e::CaloDataAnalyzer::CaloDataAnalyzer(const art::EDAnalyzer::Table& c , unpackerInstanceLabel_(config().unpackerInstanceLabel()) , verbosity_(config().verbosity()) , data_type_(config().data_type()) + , caloDAQUtil_("CaloDigiFromFragments") + , maxEventNum_(config().maxEventNum()) + , fillEmptyEvents_(config().fillEmptyEvents()) { art::ServiceHandle tfs; @@ -159,25 +163,24 @@ mu2e::CaloDataAnalyzer::CaloDataAnalyzer(const art::EDAnalyzer::Table& c tree->Branch("run", &t_run, "run/I"); tree->Branch("subrun", &t_subrun, "subrun/I"); tree->Branch("nevt", &t_nevt, "nevt/I"); - tree->Branch("dtcID", &t_dtcID, "dtcID/I"); tree->Branch( "currentDTCEventWindow", &t_currentDTCEventWindow, "currentDTCEventWindow/L"); - tree->Branch( - "currentROCEventWindow", &t_currentROCEventWindow, "currentROCEventWindow[6]/L"); tree->Branch("nhits", &t_nhits, "nhits/I"); + tree->Branch("dtcID", &t_dtcID, "dtcID[nhits]/I"); tree->Branch("boardID", &t_boardID, "boardID[nhits]/I"); tree->Branch("linkID", &t_linkID, "linkID[nhits]/I"); tree->Branch("chanID", &t_chanID, "chanID[nhits]/I"); tree->Branch("errflag", &t_errflag, "errflag[nhits]/I"); tree->Branch("fff", &t_fff, "fff[nhits]/I"); - tree->Branch("time", &t_time_tot, "time[nhits]/I"); - tree->Branch("ewhit", &t_ewhit, "ewhit[nhits]/I"); + tree->Branch("timetot", &t_time_tot, "time[nhits]/I"); + tree->Branch("ewhit", &t_ewhit, "ewhit[nhits]/L"); tree->Branch("peakpos", &t_peakpos, "peakpos[nhits]/I"); tree->Branch("peakval", &t_peakval, "peakval[nhits]/I"); tree->Branch("nofsamples", &t_nofsamples, "nofsamples[nhits]/I"); tree->Branch("firstsample", &t_firstsample, "firstsample[nhits]/I"); tree->Branch("nsamples", &t_nsamples, "nsamples/I"); tree->Branch("ADC", &t_ADC, "ADC[nsamples]/I"); + tree->Branch("ADChit", &t_ADC_hit); treeHits = tfs->make("treeHits", "Hit tree"); treeHits->Branch("run", &tH_run, "run/I"); @@ -194,14 +197,19 @@ mu2e::CaloDataAnalyzer::CaloDataAnalyzer(const art::EDAnalyzer::Table& c treeHits->Branch("chanID", &tH_chanID, "chanID/I"); treeHits->Branch("errflag", &tH_errflag, "errflag/I"); treeHits->Branch("fff", &tH_fff, "fff/I"); - treeHits->Branch("time", &tH_time_tot, "time/I"); - treeHits->Branch("ewhit", &tH_ewhit, "ewhit/I"); + treeHits->Branch("timetot", &tH_time_tot, "time/I"); + treeHits->Branch("ewhit", &tH_ewhit, "ewhit/L"); treeHits->Branch("peakpos", &tH_peakpos, "peakpos/I"); treeHits->Branch("peakval", &tH_peakval, "peakval/I"); treeHits->Branch("nofsamples", &tH_nofsamples, "nofsamples/I"); treeHits->Branch("ADC", &tH_ADC); TLOG(TLVL_DEBUG + 6) << "Reading data type " << data_type_; + + total_events = 0; + total_hits = 0; + total_goodhits = 0; + total_failedhits = 0; } void mu2e::CaloDataAnalyzer::analyze(art::Event const& event) @@ -211,6 +219,9 @@ void mu2e::CaloDataAnalyzer::analyze(art::Event const& event) // (int)eventNumber << std::endl; this_eventNumber = (int)eventNumber; + if(maxEventNum_ > 0 && this_eventNumber > maxEventNum_) + return; + t_run = event.run(); t_subrun = event.subRun(); t_nevt = event.event(); @@ -218,29 +229,45 @@ void mu2e::CaloDataAnalyzer::analyze(art::Event const& event) tH_subrun = event.subRun(); tH_nevt = event.event(); - nCaloEvents = 0; - nCaloHits = 0; + t_nhits = 0; + t_nsamples = 0; + t_ADC_hit->clear(); + + nCaloEvents = 0; + nCaloHits = 0; + event_failedhits = 0; - const auto& caloDecoderColl = *event.getValidHandle( - {unpackerModuleLabel_, unpackerInstanceLabel_}); + // const auto &caloDecoderColl = + // *event.getValidHandle({unpackerModuleLabel_, + // unpackerInstanceLabel_}); - TLOG(TLVL_DEBUG + 6) << "Iterating through " << caloDecoderColl.size() << " DTCs\n"; - for(const auto& caloDTC : caloDecoderColl) + art::InputTag caloFragmentsTag_(unpackerModuleLabel_); + auto caloDecoderColl = + event.getValidHandle>( + caloFragmentsTag_); + + TLOG(TLVL_DEBUG + 6) << "Iterating through " << caloDecoderColl->size() << " DTCs\n"; + for(auto caloDTC : *caloDecoderColl) { processCaloData(caloDTC); } + if(fillEmptyEvents_ || t_nhits > 0) + { + tree->Fill(); // Only fill if we have at least 1 hit + } + total_events++; - g_eventHits->AddPoint(this_eventNumber, nCaloHits); + g_eventHits->AddPoint(this_eventNumber, t_nhits); TLOG(TLVL_DEBUG + 6) << "[CaloDataAnalyzer::analyzer] found " << nCaloEvents - << " calo subevents in event" << (int)eventNumber; - TLOG(TLVL_DEBUG + 6) << "[CaloDataAnalyzer::analyzer] found " << nCaloHits - << " calo hits in event" << (int)eventNumber; + << " calo subevents in event " << (int)eventNumber; + TLOG(TLVL_DEBUG + 6) << "[CaloDataAnalyzer::analyzer] found " << t_nhits + << " calo hits in event " << (int)eventNumber; if(nCaloEvents == 0) { TLOG(TLVL_WARNING) - << "[CaloDataAnalyzer::analyzer] found no calo subevents in event" + << "[CaloDataAnalyzer::analyzer] found no calo subevents in event " << (int)eventNumber << "!"; } @@ -251,14 +278,15 @@ void mu2e::CaloDataAnalyzer::analyze(art::Event const& event) void mu2e::CaloDataAnalyzer::processCaloData( mu2e::CalorimeterDataDecoder const& caloDecoder) { - auto& this_subevent = caloDecoder.event_; - long int thisDTCEWT = this_subevent.GetEventWindowTag().GetEventWindowTag(true); - uint64_t dtcID = this_subevent.GetDTCID(); + auto& this_subevent = caloDecoder.event_; + long int thisDTCEWT = this_subevent.GetEventWindowTag().GetEventWindowTag(true); + t_currentDTCEventWindow = thisDTCEWT; + uint64_t dtcID = this_subevent.GetDTCID(); nCaloEvents++; // Iterate over the data blocks (ROCs) std::vector dataBlocks = this_subevent.GetDataBlocks(); - uint nROCs = dataBlocks.size(); + uint nROCs = caloDecoder.block_count(); TLOG(TLVL_DEBUG + 6) << "Iterating through " << nROCs << " data blocks (ROCs)\n"; std::vector roc_hits; for(uint iroc = 0; iroc < nROCs; iroc++) @@ -291,19 +319,11 @@ void mu2e::CaloDataAnalyzer::processCaloData( else if(data_type_ == 1) { /////// DEBUG HITS /////// auto caloHits = caloDecoder.GetCalorimeterHitTestData(iroc); - uint nHits = caloHits->size(); + // auto caloHits = caloDecoder.GetCalorimeterHitTestDataNoPointer(iroc); + uint nHits = caloHits->size(); roc_hits.push_back(nHits); - t_dtcID = dtcID; - t_currentDTCEventWindow = thisDTCEWT; - t_currentROCEventWindow[nROCs] = thisROCEWT; - t_nhits = nHits; - t_nsamples = 0; - tH_dtcID = dtcID; - tH_currentDTCEventWindow = thisDTCEWT; - tH_currentROCEventWindow = thisROCEWT; - tH_nhits = nHits; - + total_hits += nHits; for(uint ihit = 0; ihit < nHits; ihit++) { mu2e::CalorimeterDataDecoder::CalorimeterHitTestDataPacket hit = @@ -311,20 +331,28 @@ void mu2e::CaloDataAnalyzer::processCaloData( std::vector hit_waveform = caloHits->at(ihit).second; // Check that the hit is good - if(hit.BeginMarker != 0xAAA) - continue; - if(hit.LastSampleMarker == 0) - continue; - if(hit_waveform.size() == 0) - continue; - if(hit.IndexOfMaxDigitizerSample >= hit_waveform.size()) + auto errorCode = caloDAQUtil_.isHitGood(caloHits->at(ihit)); + if(errorCode) + { + event_failedhits++; + total_failedhits++; + if(verbosity_ > 0) + { + std::cout << "[CaloDataAnalyzer] BAD calo hit! DTC: " << dtcID + << ", ROC: " << iroc << ", hit number: " << ihit + << " [failure code: " << errorCode << "]" << std::endl; + caloDAQUtil_.printCaloPulse(hit); + std::cout << "[CaloDataAnalyzer] \twaveform size \t" + << hit_waveform.size() << std::endl; + } continue; - + } nCaloHits++; + total_goodhits++; // Fill hists h2_channelHits->Fill(hit.BoardID, hit.ChannelID); - g_EWTs->AddPoint(thisROCEWT, hit.InPayloadEventWindowTag); + // g_EWTs->AddPoint(thisROCEWT,hit.InPayloadEventWindowTag); h1_t0->Fill(hit.Time); h1_maxIndex->Fill(hit.IndexOfMaxDigitizerSample); h1_nSamples->Fill(hit.NumberOfSamples); @@ -335,41 +363,26 @@ void mu2e::CaloDataAnalyzer::processCaloData( } // Fill trees - t_boardID[ihit] = hit.BoardID; - t_linkID[ihit] = iroc; - t_chanID[ihit] = hit.ChannelID; - t_errflag[ihit] = hit.ErrorFlags; - t_fff[ihit] = hit.LastSampleMarker; - t_time_tot[ihit] = hit.Time; - t_ewhit[ihit] = hit.InPayloadEventWindowTag; - t_peakpos[ihit] = hit.IndexOfMaxDigitizerSample; - t_peakval[ihit] = hit_waveform[hit.IndexOfMaxDigitizerSample]; - t_nofsamples[ihit] = hit.NumberOfSamples; - t_firstsample[ihit] = t_nsamples; + t_nhits++; + t_dtcID[t_nhits - 1] = dtcID; + t_boardID[t_nhits - 1] = hit.BoardID; + t_linkID[t_nhits - 1] = iroc; + t_chanID[t_nhits - 1] = hit.ChannelID; + t_errflag[t_nhits - 1] = hit.ErrorFlags; + t_fff[t_nhits - 1] = hit.LastSampleMarker; + t_time_tot[t_nhits - 1] = hit.Time; + t_ewhit[t_nhits - 1] = hit.InPayloadEventWindowTag; + t_peakpos[t_nhits - 1] = hit.IndexOfMaxDigitizerSample; + t_peakval[t_nhits - 1] = hit_waveform[hit.IndexOfMaxDigitizerSample]; + t_nofsamples[t_nhits - 1] = hit.NumberOfSamples; + t_firstsample[t_nhits - 1] = t_nsamples; for(auto adc : hit_waveform) { t_ADC[t_nsamples] = adc; t_nsamples++; } - - tH_boardID = hit.BoardID; - tH_linkID = iroc; - tH_chanID = hit.ChannelID; - tH_errflag = hit.ErrorFlags; - tH_fff = hit.LastSampleMarker; - tH_time_tot = hit.Time; - tH_ewhit = hit.InPayloadEventWindowTag; - tH_peakpos = hit.IndexOfMaxDigitizerSample; - tH_peakval = hit_waveform[hit.IndexOfMaxDigitizerSample]; - tH_nofsamples = hit.NumberOfSamples; - tH_ADC->clear(); - for(auto adc : hit_waveform) - { - tH_ADC->push_back(adc); - } - treeHits->Fill(); + t_ADC_hit->push_back(hit_waveform); } - tree->Fill(); } else if(data_type_ == 2) { /////// COUNTERS /////// @@ -398,4 +411,15 @@ void mu2e::CaloDataAnalyzer::processCaloData( } // loop over ROCs } +void mu2e::CaloDataAnalyzer::endJob() +{ + std::cout << "\n ----- [CaloDataAnalyzer] Decoding errors summary ----- " + << std::endl; + std::cout << "Total events: " << total_events << "\n"; + std::cout << "Total hits: " << total_hits << "\n"; + std::cout << "Total good hits: " << total_goodhits << "\n"; + std::cout << "Total failed hits: " << total_failedhits << " [" + << int(100. * total_failedhits / total_hits) << "%]\n"; +} + DEFINE_ART_MODULE(mu2e::CaloDataAnalyzer) diff --git a/otsdaq-mu2e-calorimeter/ArtModules/CaloDataVerifier_module.cc b/otsdaq-mu2e-calorimeter/ArtModules/CaloDataVerifier_module.cc index 0c549d4..9a3c176 100644 --- a/otsdaq-mu2e-calorimeter/ArtModules/CaloDataVerifier_module.cc +++ b/otsdaq-mu2e-calorimeter/ArtModules/CaloDataVerifier_module.cc @@ -32,33 +32,17 @@ namespace mu2e class CaloDataVerifier : public art::EDFilter { public: - struct Config - { - fhicl::Atom verbosity{ - fhicl::Name("verbosity"), fhicl::Comment("Verbosity [0-2]"), 0}; - fhicl::Atom data_type{ - fhicl::Name("dataType"), - fhicl::Comment("Data type (0:standard, 1:debug, 2:counters)"), - 0}; - fhicl::Atom metrics_level{ - fhicl::Name("metricsLevel"), fhicl::Comment("Metrics reporting level"), 1}; - fhicl::Atom produce_calo_decoders{ - fhicl::Name("produceCaloDecoders"), - fhicl::Comment("Produce calo decoders [default: true]"), - true}; - fhicl::Atom stop_on_failure{ - fhicl::Name("stopOnFailure"), - fhicl::Comment("Throw exception if checks fail [default: false]"), - false}; - fhicl::Atom check_ewts{ - fhicl::Name("checkEWTs"), - fhicl::Comment("Check for EWT continuity across events [default: false]"), - false}; - fhicl::Atom subsystem_override{ - fhicl::Name("subsystemOverride"), - fhicl::Comment("Override calo subsystem [\"calo\", \"tracker\"]"), - "calo"}; - }; + // clang-format off + struct Config { + fhicl::Atom verbosity {fhicl::Name("verbosity" ) , fhicl::Comment("Verbosity [0-2]"), 0}; + fhicl::Atom data_type {fhicl::Name("dataType" ) , fhicl::Comment("Data type (0:standard, 1:debug, 2:counters)"), 0}; + fhicl::Atom metrics_level {fhicl::Name("metricsLevel" ) , fhicl::Comment("Metrics reporting level"), 1}; + fhicl::Atom produce_calo_decoders {fhicl::Name("produceCaloDecoders" ) , fhicl::Comment("Produce calo decoders [default: true]"), true}; + fhicl::Atom stop_on_failure {fhicl::Name("stopOnFailure" ) , fhicl::Comment("Throw exception if checks fail [default: false]"), false}; + fhicl::Atom check_ewts {fhicl::Name("checkEWTs" ) , fhicl::Comment("Check for EWT continuity across events [default: false]"), false}; + fhicl::Atom subsystem_override {fhicl::Name("subsystemOverride" ) , fhicl::Comment("Override calo subsystem [\"calo\", \"tracker\"]"), "calo"}; + }; + // clang-format on explicit CaloDataVerifier(const art::EDFilter::Table& config); @@ -611,7 +595,7 @@ void mu2e::CaloDataVerifier::processCaloData( std::cout << "========== THIS EVENT HAS SOME FAILURES ===========\n"; std::cout << "===================================================\n"; std::cout << "\n----- Dumping previous event -----\n\n"; - printEvent(*previousDTCEvent); + // printEvent(*previousDTCEvent); std::cout << "\n----- Dumping current event -----\n\n"; printEvent(dtcevent); } @@ -767,7 +751,7 @@ void mu2e::CaloDataVerifier::printEvent(const DTCLib::DTC_Event& dtcevent) std::vector subevents = dtcevent.GetSubEvents(); for(uint idtc = 0; idtc < subevents.size(); idtc++) { - std::cout << "-- DTC " << idtc << " --\n"; + std::cout << "-- DTC " << int(subevents[idtc].GetDTCID()) << " --\n"; printSubEvent(subevents[idtc]); } } @@ -828,6 +812,17 @@ void mu2e::CaloDataVerifier::printROC(const DTCLib::DTC_DataBlock& dataBlock) if(word % 16 == 15) std::cout << std::endl; } + std::cout << "--- DECODED WORDS --- " << std::endl; + uint n12bitWords = 21 * (dataBlock.GetHeader()->GetPacketCount() / 2); + mu2e::CalorimeterDataDecoder::Data12bitReader reader( + reinterpret_cast(dataBlock.GetData())); + for(size_t word = 0; word < n12bitWords; word++) + { + std::cout << std::hex << std::setw(3) << std::setfill('0') << reader[word] + << std::dec << " "; + if(word % 21 == 20) + std::cout << std::endl; + } return; } diff --git a/otsdaq-mu2e-calorimeter/ArtModules/CaloiercAnalyzer_module.cc b/otsdaq-mu2e-calorimeter/ArtModules/CaloiercAnalyzer_module.cc new file mode 100644 index 0000000..c19d470 --- /dev/null +++ b/otsdaq-mu2e-calorimeter/ArtModules/CaloiercAnalyzer_module.cc @@ -0,0 +1,494 @@ +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" + +#include "Offline/RecoDataProducts/inc/CaloDigi.hh" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include "TH1.h" +#include "TH2.h" +#include "art_root_io/TFileService.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // Libreria per manipolare il formato di output +#include +#include + +namespace mu2e +{ +class CaloiercAnalyzer : public art::EDAnalyzer +{ + public: + struct Config + { + fhicl::Atom caloDigiModuleLabel{ + fhicl::Name("caloDigiModuleLabel"), + fhicl::Comment("caloDigiModuleLabel"), + ""}; + fhicl::Atom caloDigiInstanceLabel{ + fhicl::Name("caloDigiInstanceLabel"), + fhicl::Comment("caloDigiInstanceLabel"), + ""}; + fhicl::Atom verbosity{ + fhicl::Name("verbosity"), fhicl::Comment("Verbosity [0-2]"), 0}; + fhicl::Atom splineFilename{ + fhicl::Name("splineFilename"), fhicl::Comment("splineFilename"), ""}; + fhicl::Atom uset0{ + fhicl::Name("uset0"), + fhicl::Comment("Use t0 instead of fitting with templates"), + false}; + fhicl::Atom skipAfterN{ + fhicl::Name("skipAfterN"), fhicl::Comment("Don't fit after N hits"), -1}; + }; + + explicit CaloiercAnalyzer(const art::EDAnalyzer::Table& config); + void analyze(art::Event const& event) override; + void endJob() override; + + private: + std::string caloDigiModuleLabel_; + std::string caloDigiInstanceLabel_; + int verbosity_; + std::string splineFilename_; + bool uset0_; + int skipAfterN_; + + std::map> digiMap; + std::map> timeMap; + std::map> chi2Map; + std::map> chi2rMap; + + int total_events; + int last_event; + int total_hits; + int total_unsorted; + int total_badfits; + + art::ServiceHandle tfs; + art::TFileDirectory* sameChanDir; + art::TFileDirectory* sameBoardDir; + art::TFileDirectory* diffBoardDir; + + std::map map_h1_dt_singlechan; + std::map map_g_dtevt_singlechan; + + std::map map_h1_dt_sameboard; + std::map map_g_dtevt_sameboard; + + std::map map_h1_dt_diffboard; + std::map map_g_dtevt_diffboard; + + TSpline3* spline; + TF1* f_spline; +}; +} // namespace mu2e + +mu2e::CaloiercAnalyzer::CaloiercAnalyzer(const art::EDAnalyzer::Table& config) + : art::EDAnalyzer{config} + , caloDigiModuleLabel_(config().caloDigiModuleLabel()) + , caloDigiInstanceLabel_(config().caloDigiInstanceLabel()) + , verbosity_(config().verbosity()) + , splineFilename_(config().splineFilename()) + , uset0_(config().uset0()) + , skipAfterN_(config().skipAfterN()) +{ + TFile* filetemp = TFile::Open(splineFilename_.c_str()); + if(!filetemp->IsOpen()) + { + std::cout << splineFilename_ << " NOT FOUND" << std::endl; + } + + char hname[200]; + sprintf(hname, "spline_0_1"); + if(filetemp->Get(hname)) + spline = (TSpline3*)filetemp->Get(hname); + else + { + std::cout << "No template.." << std::endl; + } + + f_spline = new TF1( + "f_spline", + [this](double* x, double* par) { + return par[0] * this->spline->Eval(x[0] - par[1]) + par[2]; + }, + 0., + 19., + 3); + f_spline->SetParNames("scale", "tpeak", "ped"); + f_spline->SetNpx(10000); + f_spline->SetRange(spline->GetXmin(), spline->GetXmax()); + + total_events = 0; + total_hits = 0; + total_unsorted = 0; + total_badfits = 0; + // TVirtualFitter::SetDefaultFitter("Minuit"); + + art::TFileDirectory t_sameChanDir = tfs->mkdir("sameChannel"); + art::TFileDirectory t_sameBoardDir = tfs->mkdir("sameBoard"); + art::TFileDirectory t_diffBoardDir = tfs->mkdir("diffBoard"); + sameChanDir = &t_sameChanDir; + sameBoardDir = &t_sameBoardDir; + diffBoardDir = &t_diffBoardDir; +} + +void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) +{ + art::EventNumber_t eventNumber = event.event(); + int this_eventNumber = (int)eventNumber; + total_events++; + last_event = this_eventNumber; + + const auto& caloDigis = + *event.getValidHandle(consumes(caloDigiModuleLabel_)); + + // Fill time-ordered map for each sipm + digiMap.clear(); + for(uint ihit = 0; ihit < caloDigis.size(); ihit++) + { + int thisID = caloDigis[ihit].SiPMID(); + float thisTime = caloDigis[ihit].t0(); + // Check if this hit is not the last, if so insert + bool sorted = true; + for(uint storedHit = 0; storedHit < digiMap[thisID].size(); storedHit++) + { + auto storedDigi = caloDigis[digiMap[thisID][storedHit]]; + if(thisTime < storedDigi.t0()) + { + std::cout << "Found unsorted hit! (in position " << storedHit << " of " + << digiMap[thisID].size() << ") "; + std::cout << "Current t0: " << thisTime + << " , last t0: " << storedDigi.t0() << "\n"; + digiMap[thisID].insert(digiMap[thisID].begin() + storedHit, ihit); + sorted = false; + total_unsorted++; + break; + } + } + if(sorted) + { // this triggers for the first one and when all is sorted + digiMap[thisID].push_back(ihit); + } + total_hits++; + } + + // Check that we filled the map properly + if(verbosity_ > 0) + { + std::cout << "\n-- LIST OF CALODIGIS READ --\n"; + for(auto pair : digiMap) + { + int thisID = pair.first; + std::cout << "SiPM ID " << thisID << " [dtc: " << thisID / 120 + << ", roc: " << (thisID % 120) / 20 << " (" << (thisID / 20) % 6 + << "), chan: " << thisID % 20 << "] (" << pair.second.size() + << " hits) : "; + for(int idx : pair.second) + { + float thisTime = caloDigis[idx].t0(); + std::cout << thisTime << " "; + } + std::cout << "\n"; + } + } + + // Now, fit the templates and fill a timeMap + timeMap.clear(); + chi2Map.clear(); + chi2rMap.clear(); + if(uset0_) + { // Use t0 instead of fitting + for(auto pair : digiMap) + { + for(int idx : pair.second) + { + timeMap[pair.first].push_back(5. * caloDigis[idx].t0()); + } + } + } + else + { // Do template fitting + // TGraphErrors *gadc = new TGraphErrors(); + TGraph* gadc = new TGraph(); + for(auto pair : digiMap) + { + for(uint ihit = 0; ihit < pair.second.size(); ihit++) + { + if(skipAfterN_ > 0 && int(ihit) >= skipAfterN_) + continue; // stop fitting if surpassed N hits + uint idx = pair.second[ihit]; // idx is the index in the original + // caloDigis vector + gadc->Set(0); + // std::cout<<"Fitting sipm "<SetPoint(gi, gi, this_waveform[gi]); + // gadc->SetPointError(gi, 0., 1.); + } + // std::cout<<" ("<SetParameters(3850, 0, 0); + // double xmin, xmax; + // f_spline->GetRange(xmin, xmax); + // std::cout<<"Now fitting in range ("<Fit("f_spline", "QRN"); // FIT! + if(fitStatus >= 0) + { + // std::cout<<"fitStatus: "<GetParameter(2)<<", time: + // "<GetParameter(1)<<", scale: + // "<GetParameter(0)<<"\n"; std::cout<<"chi2: + // "<GetChisquare()<<", ndf: "<GetNDF()<<", + // chi2/ndf: "<GetChisquare()/f_spline->GetNDF()<<"\n"; + timeMap[pair.first].push_back( + 5. * (caloDigis[idx].t0() + f_spline->GetParameter(1))); + chi2Map[pair.first].push_back(f_spline->GetChisquare()); + chi2rMap[pair.first].push_back(f_spline->GetChisquare() / + f_spline->GetNDF()); + } + else + { // if bad fit, don't fill at all + std::cout << "Bad fit status: " << fitStatus << " for sipmid " + << pair.first << " hit " << idx << "\n"; + total_badfits++; + } + } + } + delete gadc; + } + + // Print fit results + if(verbosity_ > 0) + { + std::cout << "\n-- LIST OF PULSES FITTED --\n"; + for(auto pair : timeMap) + { + std::cout << "SiPMID " << pair.first << " \t: "; + for(float this_time : pair.second) + { + std::cout << this_time << " "; + } + std::cout << "\n"; + std::cout << "Chi2 " << pair.first << " \t: "; + for(float this_chi2 : chi2Map[pair.first]) + { + std::cout << this_chi2 << " "; + } + std::cout << "\n"; + std::cout << "Chi2r " << pair.first << " \t: "; + for(float this_chi2 : chi2rMap[pair.first]) + { + std::cout << this_chi2 << " "; + } + std::cout << "\n"; + } + } + + //--------Now we can fill hists + + // Same channel + for(auto pair : timeMap) + { + int thisID = pair.first; + if(pair.second.size() > 1) + { // there must be at least 2 hits + if(map_h1_dt_singlechan.find(thisID) == map_h1_dt_singlechan.end()) + { // This hist doesn't exist yet + TString hname = Form("h1_dt_schan_%d", thisID); + TString htitle = + Form("[SAME CHANNEL] dt between hit 0 and hit 1 of sipm %d", thisID); + map_h1_dt_singlechan[thisID] = + tfs->make(hname, htitle, 600, 10009, 10011); + map_h1_dt_singlechan[thisID]->GetXaxis()->SetTitle("dt [ns]"); + TString gname = Form("g_dtevt_schan_%d", thisID); + TString gtitle = + Form("[SAME CHANNEL] dt between hit 0 and hit 1 of sipm %d", thisID); + map_g_dtevt_singlechan[thisID] = + tfs->makeAndRegister(gname, gtitle); + map_g_dtevt_singlechan[thisID]->GetXaxis()->SetTitle("Event number"); + map_g_dtevt_singlechan[thisID]->GetYaxis()->SetTitle("dt [ns]"); + map_g_dtevt_singlechan[thisID]->SetMarkerStyle(20); + } + float dt = timeMap[thisID][1] - timeMap[thisID][0]; + map_h1_dt_singlechan[thisID]->Fill(dt); + map_g_dtevt_singlechan[thisID]->AddPoint(this_eventNumber, dt); + } + } + + // Same board + int previous_DTCROC = -1; + for(auto pair : timeMap) + { + int thisID = pair.first; + int thisDTCROC = thisID / 20; + int thisDTC = thisDTCROC / 6; + int thisROC = thisDTCROC % 6; + int thisChan = thisID % 20; + if(thisDTCROC == previous_DTCROC) + continue; // Skip if we are on the same board as previous + + previous_DTCROC = thisDTCROC; + for(auto nextpair : timeMap) + { + int nextID = nextpair.first; + int nextDTCROC = nextID / 20; + int nextChan = nextID % 20; + if(thisID == nextID) + continue; + if(thisDTCROC == nextDTCROC) + { // We found a different channel in the same board + int idpair = thisID * 10000 + nextID; + if(map_h1_dt_sameboard.find(idpair) == map_h1_dt_sameboard.end()) + { // This hist doesn't exist yet + TString hname = Form("map_h1_dt_sameboard_%d_%d", nextID, thisID); + TString htitle = Form( + "[SAME BOARD] dt between chan %d and %d (DTC: %d, Board: %d, hit " + "0)", + nextChan, + thisChan, + thisDTC, + thisROC); + map_h1_dt_sameboard[idpair] = + tfs->make(hname, htitle, 5000, -50, 50); + map_h1_dt_sameboard[idpair]->GetXaxis()->SetTitle("dt [ns]"); + TString gname = Form("map_g_dtevt_sameboard_%d_%d", nextID, thisID); + TString gtitle = Form( + "[SAME BOARD] dt between chan %d and %d (DTC: %d, Board: %d, hit " + "0)", + nextChan, + thisChan, + thisDTC, + thisROC); + map_g_dtevt_sameboard[idpair] = + tfs->makeAndRegister(gname, gtitle); + map_g_dtevt_sameboard[idpair]->GetXaxis()->SetTitle("Event number"); + map_g_dtevt_sameboard[idpair]->GetYaxis()->SetTitle("dt [ns]"); + map_g_dtevt_sameboard[idpair]->SetMarkerStyle(20); + } + float dt = timeMap[nextID][0] - timeMap[thisID][0]; + if(abs(dt) > 5000) + continue; // They are not the same hit, discard + // std::cout<<"Same board, filling evt "<Fill(dt); + map_g_dtevt_sameboard[idpair]->AddPoint(this_eventNumber, dt); + } + } + } + + // Different board, same chan + // std::cout<make(hname, htitle, 20000, -200, 200); + map_h1_dt_diffboard[idpair]->GetXaxis()->SetTitle("dt [ns]"); + TString gname = Form("map_g_dtevt_diffboard_chan%d_%d_%d", + thisChan, + thisID, + refID); + TString gtitle = Form( + "[DIFFERENT BOARD] dt between (DTC: %d, Board: %d) and (DTC: " + "%d, Board: %d) [chan %d, hit 0]", + thisDTC, + thisROC, + refDTC, + refROC, + thisChan); + map_g_dtevt_diffboard[idpair] = + tfs->makeAndRegister(gname, gtitle); + map_g_dtevt_diffboard[idpair]->GetXaxis()->SetTitle( + "Event number"); + map_g_dtevt_diffboard[idpair]->GetYaxis()->SetTitle("dt [ns]"); + map_g_dtevt_diffboard[idpair]->SetMarkerStyle(20); + } + float dt = timeMap[thisID][0] - timeMap[refID][0]; + if(abs(dt) > 5000) + continue; // They are not the same hit, discard + map_h1_dt_diffboard[idpair]->Fill(dt); + map_g_dtevt_diffboard[idpair]->AddPoint(this_eventNumber, dt); + } + } + } + } +} + +void mu2e::CaloiercAnalyzer::endJob() +{ + std::cout << "\n-- JOB SUMMARY --\n"; + std::cout << "Total events: " << total_events << " (last event: " << last_event + << ")\n"; + std::cout << "Total hits: " << total_hits << "\n"; + std::cout << "Total unsorted hits found: " << total_unsorted << "\n"; + std::cout << "Total bad fits: " << total_badfits << "\n"; +} + +DEFINE_ART_MODULE(mu2e::CaloiercAnalyzer) diff --git a/test/generateCaloDigiFromFragment.fcl b/test/generateCaloDigiFromFragment.fcl new file mode 100644 index 0000000..97915de --- /dev/null +++ b/test/generateCaloDigiFromFragment.fcl @@ -0,0 +1,114 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name : FragmentToDigi + +source : { + module_type : RootInput + fileNames : @nil + maxEvents : -1 + inputCommands: [ + "keep *", + "drop *_*STM*_*_*" + ] +} + +services : { + + #@local::Services.Reco + GeometryService : {inputFile: "Offline/Mu2eG4/geom/geom_common.txt" bFieldFile: "Offline/Mu2eG4/geom/bfgeom_reco_v01.txt" simulatedDetector : {tool_type: "Mu2e"}} + ConditionsService : { conditionsfile : "Offline/ConditionsService/data/conditions_01.txt"} + GlobalConstantsService : {inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt"} + DbService : @local::DbEmpty + ProditionsService : @local::Proditions + TimeTracker: {} +} + +physics : { + + filters : { + prescaler : { + module_type : EventPrescaler + prescaleFactor : 100 + } + + caloVerifier : { + module_type : CaloDataVerifier + verbosity : 2 + dataType : 1 + metricsLevel : 2 + subsystemOverride : calo + produceCaloDecoders : true + stopOnFailure : false + checkEWTs : false + } + } + + producers : { + genFrags : { + module_type : ArtFragmentsFromDTCEvents + diagLevel : 0 + makeCaloFrag : 1 + makeTrkFrag : 0 + makeCRVFrag : 0 + makeSTMFrag : 0 + } + + CaloDigi : { + module_type : CaloDigiFromFragments + dataType : 1 + diagLevel : 0 + caloTag : "genFrags" + useOfflineID : true + } + + CaloDigiIERC : { + module_type : CaloDigiFromFragments + dataType : 1 + diagLevel : 0 + caloTag : "genFrags" + useOfflineID : false + useDTCROCID : true + } + + + } + + analyzers : { + + caloAnalyzer : { + module_type : CaloDataAnalyzer + unpackerModuleLabel : "genFrags" + unpackerInstanceLabel : "" + verbosity : 0 + dataType : 1 + } + } + + t0 : [caloVerifier] + + + t1 : [ genFrags , CaloDigiIERC ] + e1 : [ outfile, caloAnalyzer ] + + + t2 : [ genFragsSidet ] + e2 : [ caloAnalyzerSidet ] + + trigger_paths : [t1] + end_paths : [e1] +} + +outputs: { + outfile : { + module_type : RootOutput + fileName : "digis_from_frags_debug.art" + + outputCommands: [ + "drop *_*_*_*", + "keep *_*_*_*" + ] + } +} + +services.scheduler.wantSummary: true diff --git a/test/generateCaloDigiFromFragment_Sidet.fcl b/test/generateCaloDigiFromFragment_Sidet.fcl new file mode 100644 index 0000000..3dacad8 --- /dev/null +++ b/test/generateCaloDigiFromFragment_Sidet.fcl @@ -0,0 +1,96 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name : FragmentToDigi + +source : { + module_type : RootInput + fileNames : @nil + maxEvents : -1 +} + +services : { + + #@local::Services.Reco + GeometryService : {inputFile: "Offline/Mu2eG4/geom/geom_common.txt" bFieldFile: "Offline/Mu2eG4/geom/bfgeom_reco_v01.txt" simulatedDetector : {tool_type: "Mu2e"}} + #ConditionsService : { conditionsfile : "Offline/ConditionsService/data/conditions_01.txt"} + GlobalConstantsService : {inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt"} + DbService : @local::DbEmpty + ProditionsService : @local::Proditions + TimeTracker: {} +} + +physics : { + + filters : { + + caloVerifier : { + module_type : CaloDataVerifier + verbosity : 2 + dataType : 1 + metricsLevel : 2 + subsystemOverride : calo + produceCaloDecoders : true + stopOnFailure : false + checkEWTs : false + } + } + + producers : { + genFrags : { + module_type : ArtFragmentsFromDTCEvents + diagLevel : 0 + makeCaloFrag : 1 + makeTrkFrag : 0 + makeCRVFrag : 0 + makeSTMFrag : 0 + } + + CaloDigi : { + module_type : CaloDigiFromFragments + dataType : 1 + diagLevel : 0 + caloTag : "genFrags" + useOfflineID : true + } + + } + + analyzers : { + + caloAnalyzer : { + module_type : CaloDataAnalyzer + unpackerModuleLabel : "genFrags" + unpackerInstanceLabel : "" + verbosity : 0 + dataType : 1 + } + + } + + t0 : [caloVerifier] + e0 : [ ] + + t1 : [ genFrags , CaloDigi ] + e1 : [ outfile ] + + t2 : [ genFrags ] + e2 : [ caloAnalyzer ] + + trigger_paths : [t1] + end_paths : [e1] +} + +outputs: { + outfile : { + module_type : RootOutput + fileName : "digis_from_frags_debug.art" + + outputCommands: [ + "drop *_*_*_*", + "keep *_*_*_*" + ] + } +} + +services.scheduler.wantSummary: true diff --git a/test/iercAnalysis.fcl b/test/iercAnalysis.fcl new file mode 100644 index 0000000..452a123 --- /dev/null +++ b/test/iercAnalysis.fcl @@ -0,0 +1,32 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name : iercAnalysis + +source : { + module_type : RootInput + fileNames : @nil + maxEvents : -1 +} + + +physics : { + + analyzers : { + pulseIERCAnalysis : { + module_type : CaloiercAnalyzer + caloDigiModuleLabel : "CaloDigiIERC" + caloDigiInstanceLabel : "" + verbosity : 0 + splineFilename : "/home/mu2edaq/dev_calo/pgirotti_testing/srcs/otsdaq-mu2e-calorimeter/test/ierc_splines_gr4_106139.root" + } + } + + t1 : [ ] + e1 : [ pulseIERCAnalysis ] + + trigger_paths : [] + end_paths : [e1] +} + +services.scheduler.wantSummary: true