diff --git a/DAQ/CMakeLists.txt b/DAQ/CMakeLists.txt index af44480fa3..050dd78bf9 100644 --- a/DAQ/CMakeLists.txt +++ b/DAQ/CMakeLists.txt @@ -102,8 +102,8 @@ cet_build_plugin(PrefetchDAQData art::module Offline::RecoDataProducts ) -cet_build_plugin(STMWaveformDigisFromFragments art::module - REG_SOURCE src/STMWaveformDigisFromFragments_module.cc +cet_build_plugin(STMDigisFromFragments art::module + REG_SOURCE src/STMDigisFromFragments_module.cc LIBRARIES REG Offline::DAQ Offline::ProditionsService @@ -111,6 +111,25 @@ cet_build_plugin(STMWaveformDigisFromFragments art::module Offline::RecoDataProducts ) +cet_build_plugin(STMBinaryDigisFromFragments art::module + REG_SOURCE src/STMBinaryDigisFromFragments_module.cc + LIBRARIES REG + Offline::DAQ + Offline::ProditionsService + Offline::DataProducts + Offline::RecoDataProducts +) + +cet_build_plugin(STMPrintFragments art::module + REG_SOURCE src/STMPrintFragments_module.cc + LIBRARIES REG + Offline::DAQ + Offline::ProditionsService + Offline::DataProducts + Offline::RecoDataProducts + BTrk_difAlgebra +) + cet_build_plugin(StrawDigisFromArtdaqFragments art::module REG_SOURCE src/StrawDigisFromArtdaqFragments_module.cc LIBRARIES REG @@ -195,6 +214,8 @@ cet_build_plugin(MTPHitsFromDTCEvents art::module Offline::RecoDataProducts ) +install(DIRECTORY test DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/Offline/DAQ) install_source(SUBDIRS src) install_headers(USE_PROJECT_NAME SUBDIRS inc) install_fhicl(SUBDIRS fcl SUBDIRNAME Offline/DAQ/fcl) +install_fhicl(SUBDIRS test GLOB SUBDIRNAME Offline/DAQ/test) diff --git a/DAQ/src/STMBinaryDigisFromFragments_module.cc b/DAQ/src/STMBinaryDigisFromFragments_module.cc new file mode 100644 index 0000000000..d3c0bf4304 --- /dev/null +++ b/DAQ/src/STMBinaryDigisFromFragments_module.cc @@ -0,0 +1,471 @@ +// ===================================================================== +// +// STMBinaryDigisFromFragments: Binary file writing +// +// ====================================================================== + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "fhiclcpp/ParameterSet.h" + +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/RecoDataProducts/inc/STMWaveformDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" +#include "art/Framework/Principal/Handle.h" +#include "artdaq-core-mu2e/Overlays/STMFragment.hh" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace art +{ + class STMBinaryDigisFromFragments; +} + +using art::STMBinaryDigisFromFragments; + +class art::STMBinaryDigisFromFragments : public EDProducer +{ +public: + struct Config { + fhicl::Atom stmTag {fhicl::Name("stmTag"), fhicl::Comment("Input module")}; + fhicl::Atom rawFile {fhicl::Name("rawFile"), "raw.bin"}; + fhicl::Atom zsFile {fhicl::Name("zsFile"), "zs.bin"}; + fhicl::Atom phFile {fhicl::Name("phFile"), "ph.bin"}; + fhicl::Atom rawHeaderFile {fhicl::Name("rawHeaderFile"), "rawWithHeader.bin"}; + fhicl::Atom eventFile {fhicl::Name("eventFile"), "event.bin"}; + fhicl::OptionalAtom verbosityLevel{fhicl::Name("verbosityLevel"), fhicl::Comment("Verbosity level")}; + }; + + explicit STMBinaryDigisFromFragments(const art::EDProducer::Table& config); // constructor created, config via fcl + virtual ~STMBinaryDigisFromFragments(); //declares destructor + + virtual void produce(Event &) override; + void endJob() override;//Final prinout summary + +private: + art::InputTag _stmFragmentsTag; + std::ofstream _rawOut; //Files which can be acessesed throughout + std::ofstream _zsOut; + std::ofstream _phOut; + std::ofstream _rawHeaderOut; + std::ofstream _eventOut; + + //Metrics + size_t _totalEvents{0}; //Another way to initialize to zero + size_t _totalFragments{0}; + size_t _totalContainers{0}; + size_t _totalInner{0}; + size_t _totalRaw{0}; + size_t _totalZS{0}; + size_t _totalPH{0}; + + //Additional metrics + size_t _totalZeroRaw{0}; + size_t _totalZeroZS{0}; + size_t _totalZeroPH{0}; + size_t _totalEmptyRaw{0}; + size_t _totalEmptyZS{0}; + size_t _totalEmptyPH{0}; + + //fhicl varibales + int _verbosityLevel = 0; +}; // STMDigisFromFragments + +// ====================================================================== + + +STMBinaryDigisFromFragments::STMBinaryDigisFromFragments(const art::EDProducer::Table& config) + : art::EDProducer{config} + ,_stmFragmentsTag(config().stmTag()) + ,_verbosityLevel(config().verbosityLevel() ? *(config().verbosityLevel()) : 0) + +{ + //turns files into binary + _rawOut.open(config().rawFile(), std::ios::binary); + _zsOut.open(config().zsFile(), std::ios::binary); + _phOut.open(config().phFile(), std::ios::binary); + _rawHeaderOut.open(config().rawHeaderFile(), std::ios::binary); + _eventOut.open(config().eventFile(), std::ios::binary); + + //Check to make sure we are reading the file + if(!_rawOut || !_zsOut || !_phOut || !_rawHeaderOut || !_eventOut){ + throw cet::exception("FILEOPEN")<< "Failed to open one or more output files\n"; + } + +} + +STMBinaryDigisFromFragments::~STMBinaryDigisFromFragments(){ + if ( _rawOut.is_open() ) + _rawOut.close(); + if ( _zsOut.is_open() ) + _zsOut.close(); + if ( _phOut.is_open() ) + _phOut.close(); + if( _rawHeaderOut.is_open() ) + _rawHeaderOut.close(); + if( _eventOut.is_open() ) + _eventOut.close(); +}//Closing files + +// ---------------------------------------------------------------------- + +void STMBinaryDigisFromFragments::produce(Event& event) +{ + + ++_totalEvents; //Increment Total Event Counter + + //std::unique_ptr raw_waveform_digis(new mu2e::STMWaveformDigiCollection); + //std::unique_ptr zs_waveform_digis(new mu2e::STMWaveformDigiCollection); + //std::unique_ptr ph_digis(new mu2e::STMPHDigiCollection); + //std::unique_ptr raw_header_waveform_digis(new mu2e::STMWaveformDigiCollection); + //std::unique_ptr ph_waveform_digis(new mu2e::STMWaveformDigiCollection);//Original + + art::Handle STMFragmentsH; + event.getByLabel(_stmFragmentsTag, STMFragmentsH); + const auto STMFragments = STMFragmentsH.product(); + + auto writePayload = [](std::ofstream& out, // [] is a capture list that only works with these parameters + const int16_t* data, + size_t words){ + out.write(reinterpret_cast(data),//Writing to file stream in binary + words * sizeof(int16_t)); + }; + + //Event Metrics + size_t localRaw_frags{0}; + size_t localZS_frags{0}; + size_t localPH_frags{0}; + size_t zeroRaw_frags{0}; + size_t zeroZS_frags{0}; + size_t zeroPH_frags{0}; + size_t emptyRaw_frags{0}; + size_t emptyZS_frags{0}; + size_t emptyPH_frags{0}; + + size_t extractedRawWaveforms{0}; + size_t extractedZSWaveforms{0}; + size_t extractedPHDigis{0}; + //uint16_t ZSfromRaw{0}; //Keep track of ZS length from raw header + //bool readRawZSinfo{false}; + + //loop over frags + for (const auto& frag : *STMFragments) { + ++_totalFragments; //Increment Total Frag counter + + if (_verbosityLevel >=3){ std::cout <<"\nFrag_ID : " << frag.fragmentID() << "\n";} + + //Check if this is a container fragment + if (frag.type() == artdaq::Fragment::ContainerFragmentType){ + + artdaq::ContainerFragment cont_frag(frag); + ++_totalContainers; + size_t blocks = cont_frag.block_count(); + _totalInner += blocks; + + //loop over container + // i index corresponds to inner frag + for (size_t i = 0; i < cont_frag.block_count(); ++i){ + + auto inner_frag = cont_frag.at(i); + mu2e::STMFragment stm_frag(*inner_frag); + mu2e::STMWaveformDigi stm_waveform; + + //auto ptr = stm_frag.payloadBegin();//Outer variable + //auto words = stm_frag.payloadWords(); + + if (stm_frag.isRaw()) { + ++_totalRaw; //Increment job counter + ++localRaw_frags;//Increment event counter + + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true;//Assume raw frag is zero filled + + //Checks for empty frag + if (payloadWords == 0) { + if(_verbosityLevel >=3){std::cout<< "\nFound an empty frag, i = " << i <<" @Raw\n";} + ++_totalEmptyRaw;//Job counter + ++emptyRaw_frags;//Event counter + continue;//stops this loop check and goes to next frag + } + + //Check if any data points are not zero + for (size_t k =0; k < payloadWords; ++k){ + if (payloadPtr[k] != 0){ + allZeros = false; + break; + } + } + + //Check if zero filled + if (allZeros){ + if (_verbosityLevel >=3){std::cout << "\nFound a zero filled frag, i = " << i << " @Raw\n";} + ++_totalZeroRaw;//Increment counters + ++zeroRaw_frags;//Increment Zero Raw Counter + continue; + } + + //Print first 20 values + if (_verbosityLevel >=3){std::cout << "\nFound a good frag, i = " << i <<" @Raw\n";} + + if (_verbosityLevel >= 4){ + std::cout << "\nFirst 20 adcs: "; + for (size_t kk = 0; kk < std::min(payloadWords,20); ++kk){ + std::cout << payloadPtr[kk] << " ,"; + } + if(_verbosityLevel >=5){ + std::cout << "\nRaw header : Raw Length = " << stm_frag.rawLength() + <<" , ZS Length = " << stm_frag.zsLength() + << " , ZS Regions = " << stm_frag.zsRegions() + << " , i = " << i << " @Raw\n"; + } + } + //Ideally only good frags get up to here + //------full data (header + payload) + //Waveform with data creation + auto dataPtr = stm_frag.dataBegin();//Inner variable + auto dataWords = stm_frag.dataWords(); + writePayload( _rawHeaderOut , dataPtr, dataWords); + + //----- payload-only + writePayload(_rawOut, payloadPtr, payloadWords); + ++extractedRawWaveforms; + + }//End of isRaw + + else if (stm_frag.isZS()){ + ++_totalZS; //Incremenet ZS counter + ++localZS_frags; + + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true; //assumes all adcs are zero + + //Check if payload is empty + if (payloadWords == 0) { + if (_verbosityLevel >=3){std::cout << "\nFound an empty frag, i = " << i << " @ZS\n";} + ++_totalEmptyZS; + ++emptyZS_frags; //Increment empty ZS counter + continue; + } + //Check if any adc are non-zero + for (size_t k = 0 ; k < payloadWords ; ++k){ + if(payloadPtr[k] != 0 ){ + allZeros = false; + break; + } + } + //Check if zero filled + if (allZeros) { + if (_verbosityLevel >=3){std::cout << "\nFound a zero filled frag, i = "<< i << " @ZS\n";} + ++_totalZeroZS; + ++zeroZS_frags; + continue; + } + //Print first 20 payload adcs + if (_verbosityLevel >=3){std::cout << "\nFound a good frag, i = " << i << " @ZS\n";} + + if (_verbosityLevel >=4){ + std::cout << "\nFirst 20 adcs: "; + for (size_t kk = 0 ; kk < std::min(payloadWords,20) ; ++kk){ + std::cout << payloadPtr[kk] << " , "; + } + std::cout << "i = " << i << " @ZS\n"; + } + + //Defintions for payload references + auto dataPtr = stm_frag.dataBegin(); + auto dataWords = stm_frag.dataWords(); + auto dataEnd = dataPtr + dataWords; + size_t seg = 0; + size_t totalLen = 0; + uint16_t lastZSindex = 0; //keeps track of last recorded index from header -> with respect to what? + uint16_t lastLen = 0; //keeps track of last recorded length from header + + writePayload(_zsOut, payloadPtr, payloadWords); + + while (dataPtr + 2 <= dataEnd){ + uint16_t current_zs_location = static_cast(dataPtr[0]); + uint16_t current_zs_size = static_cast(dataPtr[1]); + auto adc = dataPtr + 2; + if (adc + current_zs_size > dataEnd){ break;} + + ++extractedZSWaveforms; + + if (_verbosityLevel >=6){ + //A print check per segment + std::cout << "Region = " << seg << " , zs_index = " << current_zs_location << " , zs_size = " << current_zs_size + << "\n" ; + } + + //Update variables + lastZSindex = current_zs_location; + lastLen = current_zs_size; + totalLen += lastLen; + ++seg; + dataPtr = adc + current_zs_size; + } + + if (_verbosityLevel >=5){ + //Summary + std::cout << "ZS Regions = " << seg + << " , lastZSindex = " << lastZSindex + << " , lastZSLen = " << lastLen + << " , ZS total length = " << totalLen + << " , i = " << i << " @ZS\n"; + } + + //Throw out if ZSLengthfromRaw != totalLen + //if (readRawZSinfo && ZSfromRaw != totalLen){ + //throw cet::exception("Mismatch") + // << "\n=== ZS Length mismatch ===\n" + // << "ZS length from Raw header " << ZSfromRaw << "\n" + // << "ZS length calculated from file : " << totalLen << "\n" + // << "Found at inner frag i : " << i <<"\n"; + // } + + }//End of isZS + + else if (stm_frag.isPH()){ + ++_totalPH; + ++localPH_frags; + //Check if zero filled + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true; + + if (payloadWords ==0){ + if (_verbosityLevel >= 3){std::cout << "\nFound an empty frag, i = " << i<< " @PH\n";} + ++_totalEmptyPH; + ++emptyPH_frags; + continue; + } + + //Check if zero filled + for (size_t k = 0; k< payloadWords; ++k) { + if (payloadPtr[k] !=0){ + allZeros = false; + break; + } + } + //count if zero filled + if (allZeros){ + if (_verbosityLevel >=3){ std::cout<< "\nFound a zero filled frag, i = " << i<< " @PH\n";} + ++_totalZeroPH; + ++zeroPH_frags; + continue; + } + + if (_verbosityLevel >= 3){std::cout << "\nFound a good frag, i = " << i <<" @PH\n";} + + size_t digiWords = stm_frag.payloadWords(); + auto const* digiPtr = stm_frag.payloadBegin(); + + writePayload(_phOut,digiPtr, digiWords); + extractedPHDigis += digiWords; + + }//End of isPH and is checks + + //---Combined stream write w. order preserved ---- + { + const int16_t* cptr = nullptr; + size_t cwords = 0; + + if (stm_frag.isRaw()){ + cptr = stm_frag.dataBegin(); + cwords = stm_frag.dataWords(); + } + else if (stm_frag.isZS()){ + cptr = stm_frag.payloadBegin(); + cwords = stm_frag.payloadWords(); + } + else if (stm_frag.isPH()){ + cptr = stm_frag.payloadBegin(); + cwords = stm_frag.payloadWords(); + } + + if (cptr != nullptr && cwords >0) { + _eventOut.write(reinterpret_cast(cptr), + cwords * sizeof(int16_t) ); + } + } + + }//End Container loop + } else { + //fallback (non-container case) + mu2e::STMFragment stm_frag(frag); + auto ptr = stm_frag.payloadBegin(); + auto words = stm_frag.payloadWords(); + + if (stm_frag.isRaw()) { + writePayload(_rawOut, ptr, words); + } + } + } //End of frags loop---George suggestions + if (_verbosityLevel >= 2 ){ + //Event Summary -> tells us what happens per event + std::cout << "\n========== STM EVENT SUMMARY - (Binary Module) ==========\n"; + std::cout << "Extracted Raw waveforms : " << extractedRawWaveforms <<"\n"; + std::cout << "Extracted ZS waveforms : " << extractedZSWaveforms <<"\n"; + std::cout << "Extracted PH digis : " << extractedPHDigis <<"\n"; + + std::cout << "\n--- Frags Read ---\n"; + std::cout << "Raw frags : " << localRaw_frags << "\n"; + std::cout << "ZS frags : " << localZS_frags << "\n"; + std::cout << "PH frags : " << localPH_frags << "\n"; + + std::cout << "\n--- Filter results ---\n"; + std::cout << "Zero Raw frags : " << zeroRaw_frags << "\n"; + std::cout << "Zero ZS frags : " << zeroZS_frags << "\n"; + std::cout << "Zero PH frags : " << zeroPH_frags << "\n"; + + std::cout << "Empty Raw frags : " << emptyRaw_frags <<"\n"; + std::cout << "Empty ZS frags : " << emptyZS_frags << "\n"; + std::cout << "Empty PH frags : " << emptyPH_frags << "\n"; + + std::cout << "=================================\n"; + } + +} // produce() + +// ====================================================================== + +void STMBinaryDigisFromFragments::endJob() { + if (_verbosityLevel >= 1){ + //Tells us what happened at the very end + std::cout << "\n========== STM JOB SUMMARY - (Binary Module) ==========\n"; + + std::cout << "Total events : " << _totalEvents << "\n"; + std::cout << "Total fragments : " << _totalFragments << "\n"; + std::cout << "Container frags : " << _totalContainers << "\n"; + std::cout << "Inner fragments : " << _totalInner << "\n"; + + std::cout << "\n--- Data types read ---\n"; + std::cout << "RAW : " << _totalRaw << "\n"; + std::cout << "ZS : " << _totalZS << "\n"; + std::cout << "PH : " << _totalPH << "\n"; + + std::cout << "\n--- Data types filtered ---\n"; + std::cout << "Zero RAW frags : " << _totalZeroRaw << "\n"; + std::cout << "Zero ZS frags : " << _totalZeroZS << "\n"; + std::cout << "Zero PH frags : " << _totalZeroPH << "\n"; + std::cout << "Empty Raw frags : " << _totalEmptyRaw << "\n"; + std::cout << "Empty ZS frags : " << _totalEmptyZS << "\n"; + std::cout << "Empty PH frags : " << _totalEmptyPH << "\n"; + + std::cout << "=================================\n"; + + } +} + +DEFINE_ART_MODULE(STMBinaryDigisFromFragments) diff --git a/DAQ/src/STMDigisFromFragments_module.cc b/DAQ/src/STMDigisFromFragments_module.cc new file mode 100644 index 0000000000..56be40896b --- /dev/null +++ b/DAQ/src/STMDigisFromFragments_module.cc @@ -0,0 +1,683 @@ +// ===================================================================== +// +// STMDigisFromFragments: create all types of STMDigis from STMFragments +// +// ====================================================================== + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "fhiclcpp/ParameterSet.h" + +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/RecoDataProducts/inc/STMWaveformDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" +#include "art/Framework/Principal/Handle.h" +#include "artdaq-core-mu2e/Overlays/STMFragment.hh" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace art +{ + class STMDigisFromFragments; +} + +using art::STMDigisFromFragments; + +class art::STMDigisFromFragments : public EDProducer +{ +public: + struct Config { + fhicl::Atom stmTag {fhicl::Name("stmTag"), fhicl::Comment("Input module")}; + fhicl::OptionalAtom verbosityLevel{fhicl::Name("verbosityLevel"), fhicl::Comment("Verbosity level")}; + fhicl::Atom saveRawWithHeaderWaveform_HPGe {fhicl::Name("saveRawWithHeaderWaveform_HPGe"), false}; + fhicl::Atom saveRawWaveform_HPGe {fhicl::Name("saveRawWaveform_HPGe"), false}; + fhicl::Atom saveZSWaveform_HPGe{fhicl::Name("saveZSWaveform_HPGe"), false}; + fhicl::Atom saveRawWithHeaderWaveform_LaBr{fhicl::Name("saveRawWithHeaderWaveform_LaBr"), false}; + fhicl::Atom saveRawWaveform_LaBr{fhicl::Name("saveRawWaveform_LaBr"), false}; + fhicl::Atom saveZSWaveform_LaBr{fhicl::Name("saveZSWaveform_LaBr"), false}; + + }; + + explicit STMDigisFromFragments(const art::EDProducer::Table& config); // constructor created, config via fcl + virtual void produce(Event &) override; + void endJob() override;//Final prinout summary + +private: + art::InputTag _stmFragmentsTag; + + //Metrics + size_t _totalEvents{0}; + size_t _totalFragments{0}; + size_t _totalContainers{0}; + size_t _totalContainersHPGe{0}; + size_t _totalContainersLaBr{0}; + size_t _totalInner{0}; + size_t _totalRaw{0}; + size_t _totalZS{0}; + size_t _totalPH{0}; + size_t _totalZeroRaw{0}; + size_t _totalZeroZS{0}; + size_t _totalZeroPH{0}; + size_t _totalEmptyRaw{0}; + size_t _totalEmptyZS{0}; + size_t _totalEmptyPH{0}; + size_t _unreadInnerFrags{0}; + size_t _totalEventsWithHPGeLaBr{0}; + size_t _totalEventsWithOnlyHPGe{0}; + size_t _totalEventsWithOnlyLaBr{0}; + size_t _totalEventsWithNone{0}; + + //Additional metrics + //HPGe + size_t _totalRawHPGe{0}; + size_t _totalZSHPGe{0}; + size_t _totalPHHPGe{0}; + size_t _totalGoodRawHPGe{0}; + size_t _totalGoodZSHPGe{0}; + size_t _totalGoodPHHPGe{0}; + size_t _totalZeroRawHPGe{0}; + size_t _totalZeroZSHPGe{0}; + size_t _totalZeroPHHPGe{0}; + size_t _totalEmptyRawHPGe{0}; + size_t _totalEmptyZSHPGe{0}; + size_t _totalEmptyPHHPGe{0}; + + //LaBr + size_t _totalRawLaBr{0}; + size_t _totalZSLaBr{0}; + size_t _totalPHLaBr{0}; + size_t _totalGoodRawLaBr{0}; + size_t _totalGoodZSLaBr{0}; + size_t _totalGoodPHLaBr{0}; + size_t _totalZeroRawLaBr{0}; + size_t _totalZeroZSLaBr{0}; + size_t _totalZeroPHLaBr{0}; + size_t _totalEmptyRawLaBr{0}; + size_t _totalEmptyZSLaBr{0}; + size_t _totalEmptyPHLaBr{0}; + + //fhicl varibales + int _verbosityLevel = 0; + bool _saveRawWithHeaderWaveform_HPGe{false}; + bool _saveRawWaveform_HPGe{false}; + bool _saveZSWaveform_HPGe{false}; + bool _saveRawWithHeaderWaveform_LaBr{false}; + bool _saveRawWaveform_LaBr{false}; + bool _saveZSWaveform_LaBr{false}; + + +}; // STMDigisFromFragments + +// ====================================================================== + + +STMDigisFromFragments::STMDigisFromFragments(const art::EDProducer::Table& config) + : art::EDProducer{config} + ,_stmFragmentsTag(config().stmTag()) + ,_saveRawWithHeaderWaveform_HPGe(config().saveRawWithHeaderWaveform_HPGe()) + ,_saveRawWaveform_HPGe(config().saveRawWaveform_HPGe()) + ,_saveZSWaveform_HPGe(config().saveZSWaveform_HPGe()) + ,_saveRawWithHeaderWaveform_LaBr(config().saveRawWithHeaderWaveform_LaBr()) + ,_saveRawWaveform_LaBr(config().saveRawWaveform_LaBr()) + ,_saveZSWaveform_LaBr(config().saveZSWaveform_LaBr()) + ,_verbosityLevel(config().verbosityLevel() ? *(config().verbosityLevel()) : 0) + +{ + // Set the size of vectors for HPGe + if (_saveRawWithHeaderWaveform_HPGe){produces("rawWithHeaderHPGe");} + if (_saveRawWaveform_HPGe){produces("rawHPGe");}//Waveforms + if (_saveZSWaveform_HPGe){produces("zsHPGe");} + produces("phHPGe"); // digi series + + + //Set the size of vectors for LaBr + if (_saveRawWithHeaderWaveform_LaBr){produces("rawWithHeaderLaBr");} + if (_saveRawWaveform_LaBr){produces("rawLaBr");} + if (_saveZSWaveform_LaBr){produces("zsLaBr");} + produces("phLaBr"); // digi series + +} + +// ---------------------------------------------------------------------- + +void STMDigisFromFragments::produce(Event& event) +{ + + ++_totalEvents; //Increment Total Event Counter + + //HPGe + std::unique_ptr raw_HPGe_waveform_digis(new mu2e::STMWaveformDigiCollection); + std::unique_ptr zs_HPGe_waveform_digis(new mu2e::STMWaveformDigiCollection); + std::unique_ptr ph_HPGe_digis(new mu2e::STMPHDigiCollection); + std::unique_ptr raw_HPGe_header_waveform_digis(new mu2e::STMWaveformDigiCollection); + + //LaBr + std::unique_ptr raw_LaBr_waveform_digis(new mu2e::STMWaveformDigiCollection); + std::unique_ptr zs_LaBr_waveform_digis(new mu2e::STMWaveformDigiCollection); + std::unique_ptr ph_LaBr_digis(new mu2e::STMPHDigiCollection); + std::unique_ptr raw_LaBr_header_waveform_digis(new mu2e::STMWaveformDigiCollection); + + art::Handle STMFragmentsH; + event.getByLabel(_stmFragmentsTag, STMFragmentsH); + const auto& STMFragments = STMFragmentsH.product(); + + //Event Metrics + + //HPGe + size_t localRawHPGe_frags{0}; + size_t localZSHPGe_frags{0}; + size_t localPHHPGe_frags{0}; + size_t zeroRawHPGe_frags{0}; + size_t zeroZSHPGe_frags{0}; + size_t zeroPHHPGe_frags{0}; + size_t emptyRawHPGe_frags{0}; + size_t emptyZSHPGe_frags{0}; + size_t emptyPHHPGe_frags{0}; + + //LaBr + size_t localRawLaBr_frags{0}; + size_t localZSLaBr_frags{0}; + size_t localPHLaBr_frags{0}; + size_t zeroRawLaBr_frags{0}; + size_t zeroZSLaBr_frags{0}; + size_t zeroPHLaBr_frags{0}; + size_t emptyRawLaBr_frags{0}; + size_t emptyZSLaBr_frags{0}; + size_t emptyPHLaBr_frags{0}; + + //Additional + size_t unread_InnerFrags{0}; + uint16_t outerFragID{0}; + bool eventHasHPGe{false}; + bool eventHasLaBr{false}; + + //HPGe + uint16_t expectedZSLengthHPGe{0}; + uint16_t expectedZSRegionsHPGe{0}; + bool readZSinfoFromRawHeaderHPGe{false}; + + //LaBr + uint16_t expectedZSLengthLaBr{0}; + uint16_t expectedZSRegionsLaBr{0}; + bool readZSinfoFromRawHeaderLaBr{false}; + + //loop over outer frags + for (const auto& frag : *STMFragments) { + ++_totalFragments; //Increment Total Frag counter + + outerFragID = frag.fragmentID(); + if (_verbosityLevel >= 3){std::cout << "\nFrag_id : " << outerFragID << "\n";} + + + //Check if this is a container fragment + if (frag.type() == artdaq::Fragment::ContainerFragmentType){ + + mu2e::STMFragment container_frag(frag); + if ( container_frag.isHPGeContainer()){ ++_totalContainersHPGe; eventHasHPGe = true; } + else if ( container_frag.isLaBrContainer() ) { ++_totalContainersLaBr; eventHasLaBr = true; } + else if (_verbosityLevel >= 1){ + std::cout << "Encounter an unknown STM Container frag ID\n" + << "Frag ID : " << frag.fragmentID() <<"\n" + << "Event : " << _totalEvents <<"\n"; + } + + artdaq::ContainerFragment cont_frag(frag); + ++_totalContainers; + size_t blocks = cont_frag.block_count(); + _totalInner += blocks; + + + //loop over container where i corresponds to inner frag + for (size_t i = 0; i < cont_frag.block_count(); ++i){ + + auto inner_frag = cont_frag.at(i); + mu2e::STMFragment stm_frag(*inner_frag); + mu2e::STMWaveformDigi stm_waveform; + + if ( stm_frag.isRaw() ) { + + //Job Counter + ++_totalRaw; + + //Conditional Job and Event Counter + if( stm_frag.isHPGe() ){ + ++_totalRawHPGe ; + ++localRawHPGe_frags; + + } else if ( stm_frag.isLaBr() ){ + ++_totalRawLaBr ; + ++localRawLaBr_frags; + } + + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true;//Assume raw frag is zero filled + + //Checks for empty frag + if (payloadWords == 0) { + if(_verbosityLevel >=3){std::cout<< "\nFound an empty frag, i = " << i <<" @Raw\n";} + ++_totalEmptyRaw;//Job counter + //Eventcounter + if( stm_frag.isHPGe() ){ ++_totalEmptyRawHPGe; ++emptyRawHPGe_frags; } else if (stm_frag.isLaBr()){ ++_totalEmptyRawLaBr; ++emptyRawLaBr_frags; } + continue;//stops this loop check and goes to next frag + } + + //Check if any data points are not zero + for (size_t k =0; k < payloadWords; ++k){ + if (payloadPtr[k] != 0){ + allZeros = false; + break; + } + } + + //Check if zero filled + if (allZeros){ + if (_verbosityLevel >=3){std::cout << "\nFound a zero filled frag, i = " << i << " @Raw\n";} + ++_totalZeroRaw;//Increment counters + //Eventcounter + if( stm_frag.isHPGe() ){++_totalZeroRawHPGe; ++zeroRawHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalZeroRawLaBr; ++zeroRawLaBr_frags; } + continue; + } + + //Print first 20 values + if (_verbosityLevel >=3){std::cout << "\nFound a good frag, i = " << i <<" @Raw\n";} + + if (_verbosityLevel >= 4){ + std::cout << "\nFirst 20 adcs: "; + for (size_t kk = 0; kk < std::min(payloadWords,20); ++kk){ + std::cout << payloadPtr[kk] << " ,"; + } + if(_verbosityLevel >=5){ + std::cout << "\nRaw header : Raw Length = " << stm_frag.rawLength() + <<" , ZS Length = " << stm_frag.zsLength() + << " , ZS Regions = " << stm_frag.zsRegions() + << " , i = " << i << " @Raw\n"; + } + } + + //Ideally only good frags get up to here + if ( stm_frag.isHPGe() ){ + expectedZSRegionsHPGe = stm_frag.zsRegions(); + expectedZSLengthHPGe = stm_frag.zsLength(); + readZSinfoFromRawHeaderHPGe = true; + + ++ _totalGoodRawHPGe; + + if (_saveRawWithHeaderWaveform_HPGe){ + auto dataPtr = stm_frag.dataBegin();//memory check for emplace_back + auto dataWords = stm_frag.dataWords(); + stm_waveform.set_data(dataWords, dataPtr); + raw_HPGe_header_waveform_digis->emplace_back(stm_waveform); + } + if(_saveRawWaveform_HPGe) { + stm_waveform.set_data(payloadWords, payloadPtr); + raw_HPGe_waveform_digis->emplace_back(stm_waveform); + } + + } else if ( stm_frag.isLaBr() ){ + expectedZSRegionsLaBr = stm_frag.zsRegions(); + expectedZSLengthLaBr = stm_frag.zsLength(); + readZSinfoFromRawHeaderLaBr = true; + + ++ _totalGoodRawLaBr; + + if (_saveRawWithHeaderWaveform_LaBr){ + auto dataPtr = stm_frag.dataBegin(); + auto dataWords = stm_frag.dataWords(); + stm_waveform.set_data(dataWords, dataPtr); + raw_LaBr_header_waveform_digis->emplace_back(stm_waveform); + } + if(_saveRawWaveform_LaBr){ + stm_waveform.set_data(payloadWords, payloadPtr); + raw_LaBr_waveform_digis->emplace_back(stm_waveform); + } + + } + + }//End of isRaw + + + else if (stm_frag.isZS()){ + ++_totalZS; //Incremenet ZS counter + if ( stm_frag.isHPGe() ){ ++_totalZSHPGe; ++localZSHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalZSLaBr; ++localZSLaBr_frags; } + + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true; //assumes all adcs are zero + + //Check if payload is empty + if ( payloadWords == 0) { + if (_verbosityLevel >=3){std::cout << "\nFound an empty frag, i = " << i << " @ZS\n";} + ++_totalEmptyZS; + if ( stm_frag.isHPGe() ){ ++_totalEmptyZSHPGe; ++emptyZSHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalEmptyZSLaBr; ++emptyZSLaBr_frags; } + continue; + } + //Check if any adc are non-zero + for (size_t k = 0 ; k < payloadWords ; ++k){ + if(payloadPtr[k] != 0 ){ + allZeros = false; + break; + } + } + //Check if zero filled + if (allZeros) { + if (_verbosityLevel >=3){std::cout << "\nFound a zero filled frag, i = "<< i << " @ZS\n";} + ++_totalZeroZS; + if ( stm_frag.isHPGe() ){ ++_totalZeroZSHPGe; ++zeroZSHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalZeroZSLaBr; ++zeroZSLaBr_frags; } + continue; + } + //Print first 20 payload adcs + if (_verbosityLevel >=3){std::cout << "\nFound a good frag, i = " << i << " @ZS\n";} + + if (_verbosityLevel >=4){ + std::cout << "\nFirst 20 adcs: "; + for (size_t kk = 0 ; kk < std::min(payloadWords,20) ; ++kk){ + std::cout << payloadPtr[kk] << " , "; + } + std::cout << "i = " << i << " @ZS\n"; + } + + //Defintions for payload references + auto dataPtr = stm_frag.dataBegin(); + auto dataWords = stm_frag.dataWords(); + auto dataEnd = dataPtr + dataWords; + size_t seg = 0; + size_t totalLen = 0; + uint16_t lastZSindex = 0; //keeps track of last recorded index from header -> with respect to what? + uint16_t lastLen = 0; //keeps track of last recoded length from header + + + //Try the ? and : + bool readZSinfoFromRawHeader = stm_frag.isHPGe() ? readZSinfoFromRawHeaderHPGe : readZSinfoFromRawHeaderLaBr; + uint16_t expectedZSRegions = stm_frag.isHPGe() ? expectedZSRegionsHPGe : expectedZSRegionsLaBr; + uint16_t expectedZSLength = stm_frag.isHPGe() ? expectedZSLengthHPGe : expectedZSLengthLaBr; + + if(_verbosityLevel >= 6){std::cout << "dataWords : " << dataWords + << " dataWords%4 : " << dataWords%4 + <<" @ZS" << "\n"; } + + if ( stm_frag.isHPGe()){ ++_totalGoodZSHPGe;} else if(stm_frag.isLaBr()){++ _totalGoodZSLaBr;} + + while (dataPtr + 2 <= dataEnd){ + if ( readZSinfoFromRawHeader && seg >= expectedZSRegions ) break; + uint16_t current_zs_location = static_cast(dataPtr[0]); + uint16_t current_zs_size = static_cast(dataPtr[1]); + auto adc = dataPtr + 2; + if (adc + current_zs_size > dataEnd) + break; + + uint32_t trigTimeOffset = current_zs_location; + std::vector segADCS(adc, adc + current_zs_size); //1D array, contains adcs to this_zs_Size - 1 + mu2e::STMWaveformDigi stm_waveform(trigTimeOffset, segADCS); //New constructore use + //emplacing + if ( stm_frag.isHPGe() && _saveZSWaveform_HPGe){ + std::vector segADCS(adc, adc + current_zs_size); //1D array, contains adcs to this_zs_size - 1 + zs_HPGe_waveform_digis->emplace_back(trigTimeOffset,segADCS); + } + else if ( stm_frag.isLaBr() && _saveZSWaveform_LaBr){ + std::vector segADCS(adc, adc + current_zs_size); + zs_LaBr_waveform_digis->emplace_back(trigTimeOffset,segADCS); + } + + if (_verbosityLevel >=6){ + + //A print check per segment + std::cout << "Region = " << seg << " , zs_index = " << current_zs_location << " , zs_size = " << current_zs_size + << " , trigTimeOffset = " << trigTimeOffset << "\n" ; + } + + //Update variables + lastZSindex = current_zs_location; + lastLen = current_zs_size; + totalLen += lastLen; + ++seg; + dataPtr = adc + current_zs_size; + } + + if (_verbosityLevel >=5){ + //Summary + std::cout << "ZS Regions = " << seg + << " , lastZSindex = " << lastZSindex + << " , lastZSLen = " << lastLen + << " , ZS total length = " << totalLen + << " , i = " << i << " @ZS\n"; + } + + //Throw out if ZSLengthfromRaw != totalLen + if (readZSinfoFromRawHeader && expectedZSLength != totalLen){ + throw cet::exception("STM_UNPACKING") + << "\n=== ZS Length mismatch ===\n" + << "ZS length from Raw header : " << expectedZSLength << "\n" + << "ZS length calculated from file : " << totalLen << "\n" + << "Found at inner frag i : " << i << "\n" + << "Encountered at event : " << _totalEvents << "\n" ; + } + + //reset variables after ZS + if (stm_frag.isHPGe()){ + readZSinfoFromRawHeaderHPGe = false; + expectedZSLengthHPGe = 0; + expectedZSRegionsHPGe = 0; + + } else if (stm_frag.isLaBr()){ + readZSinfoFromRawHeaderLaBr = false; + expectedZSLengthLaBr = 0; + expectedZSRegionsLaBr = 0; + } + + }//End of isZS + + else if ( stm_frag.isPH() ){ + + ++_totalPH; + if ( stm_frag.isHPGe() ){ ++_totalPHHPGe; ++localPHHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalPHLaBr; ++localPHLaBr_frags; } + + //Check if zero filled + auto payloadPtr = stm_frag.payloadBegin(); + auto payloadWords = stm_frag.payloadWords(); + bool allZeros = true; + + if (payloadWords ==0){ + if (_verbosityLevel >= 3){std::cout << "\nFound an empty frag, i = " << i<< " @PH\n";} + ++_totalEmptyPH; + if ( stm_frag.isHPGe() ){ ++_totalEmptyPHHPGe; ++emptyPHHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalEmptyPHLaBr; ++emptyPHLaBr_frags; } + continue; + } + + //Check if zero filled + for (size_t k = 0; k< payloadWords; ++k) { + if (payloadPtr[k] !=0){ + allZeros = false; + break; + } + } + //count if zero filled + if (allZeros){ + if (_verbosityLevel >=3){ std::cout<< "\nFound a zero filled frag, i = " << i<< " @PH\n";} + ++_totalZeroPH; + if ( stm_frag.isHPGe() ){ ++_totalZeroPHHPGe; ++zeroPHHPGe_frags; } else if ( stm_frag.isLaBr() ){ ++_totalZeroPHLaBr; ++zeroPHLaBr_frags; } + continue; + } + + if (_verbosityLevel >= 3){std::cout << "\nFound a good frag, i = " << i <<" @PH\n";} + + if ( stm_frag.isHPGe() ){ ++ _totalGoodPHHPGe; } else if ( stm_frag.isLaBr() ) { ++ _totalGoodPHLaBr; } + + size_t digiWords = stm_frag.payloadWords(); + auto const* digiPtr = stm_frag.payloadBegin(); + + for (size_t i_PH = 0; i_PH < digiWords ; ++i_PH){ + int16_t PH = digiPtr[i_PH]; + mu2e::STMPHDigi PH_digi(0, PH); + + if( stm_frag.isHPGe() ){ + ph_HPGe_digis->emplace_back(PH_digi); } + else if ( stm_frag.isLaBr() ){ + ph_LaBr_digis->emplace_back(PH_digi); } + + } + + }//End of isPH and is checks + else{ + ++ _unreadInnerFrags; //For Job summary + ++ unread_InnerFrags; //For event summary + if(_verbosityLevel >=3){ + std::cout << "Unread Inner fragment" << "\n" + << "Inner frag i : " << i << "\n" + << "Frag ID : " << inner_frag->fragmentID() << "\n" + << "Event : " <<_totalEvents << "\n"; + } + } + } + + } else { + //fallback (non-container case) + if (_verbosityLevel >=1) { + std::cout << "Found non-container STM Fragment " <<"\n" + << "Fragment ID : " << frag.fragmentID() << "\n" + << "Event : " << _totalEvents << "\n"; + } + } + } //End of frags loop + + if (eventHasHPGe && eventHasLaBr){ ++_totalEventsWithHPGeLaBr; } + else if (eventHasHPGe){ ++_totalEventsWithOnlyHPGe; } + else if (eventHasLaBr){ ++_totalEventsWithOnlyLaBr; } + else { ++_totalEventsWithNone; } + + if (_verbosityLevel >= 2 ){ + //Event Summary -> tells us what happens per event + std::cout << "\n========== STM EVENT SUMMARY - (Unpacking Module) ==========\n"; + std::cout << "\n--- Products extracted ---\n"; + std::cout << "Extracted Raw HPGe waveforms : " << raw_HPGe_waveform_digis->size() <<"\n"; + std::cout << "Extracted ZS HPGe waveforms : " << zs_HPGe_waveform_digis->size() <<"\n"; + std::cout << "Extracted PH HPGe digis : " << ph_HPGe_digis->size() <<"\n"; + std::cout << "\n"; + std::cout << "Extracted Raw LaBr waveforms : " << raw_LaBr_waveform_digis->size() <<"\n"; + std::cout << "Extracted ZS LaBr waveforms : " << zs_LaBr_waveform_digis->size() <<"\n"; + std::cout << "Extracted PH LaBr digis : " << ph_LaBr_digis->size() <<"\n"; + + + std::cout << "\n--- Frags Seen ---\n"; + std::cout << "Raw HPGE frags seen : " << localRawHPGe_frags << "\n"; + std::cout << "ZS HPGe frags seen : " << localZSHPGe_frags << "\n"; + std::cout << "PH HPGe frags seen : " << localPHHPGe_frags << "\n"; + std::cout << "\n"; + std::cout << "Raw LaBr frags seen : " << localRawLaBr_frags << "\n"; + std::cout << "ZS LaBr frags seen : " << localZSLaBr_frags << "\n"; + std::cout << "PH LaBr frags seen : " << localPHLaBr_frags << "\n"; + std::cout << "\nUnread frags : " << unread_InnerFrags << "\n"; + + + std::cout << "\n--- Filter results ---\n"; + std::cout << "Zero Raw HPGe frags : " << zeroRawHPGe_frags << "\n"; + std::cout << "Zero ZS HPGe frags : " << zeroZSHPGe_frags << "\n"; + std::cout << "Zero PH HPGe frags : " << zeroPHHPGe_frags << "\n"; + std::cout << "\n"; + std::cout << "Zero Raw LaBr frags : " << zeroRawLaBr_frags << "\n"; + std::cout << "Zero ZS LaBr frags : " << zeroZSLaBr_frags << "\n"; + std::cout << "Zero PH LaBr frags : " << zeroPHLaBr_frags << "\n"; + std::cout << "\n"; + + std::cout << "Empty Raw HPGe frags : " << emptyRawHPGe_frags <<"\n"; + std::cout << "Empty ZS HPGe frags : " << emptyZSHPGe_frags << "\n"; + std::cout << "Empty PH HPGe frags : " << emptyPHHPGe_frags << "\n"; + std::cout << "\n"; + std::cout << "Empty Raw LaBr frags : " << emptyRawLaBr_frags <<"\n"; + std::cout << "Empty ZS LaBr frags : " << emptyZSLaBr_frags << "\n"; + std::cout << "Empty PH LaBr frags : " << emptyPHLaBr_frags << "\n"; + + std::cout << "=================================\n"; + } + + //Final move + //HPGe + if (_saveRawWithHeaderWaveform_HPGe){ event.put(std::move(raw_HPGe_header_waveform_digis), "rawWithHeaderHPGe"); } + if (_saveRawWaveform_HPGe){event.put(std::move(raw_HPGe_waveform_digis), "rawHPGe");} + if (_saveZSWaveform_HPGe){event.put(std::move(zs_HPGe_waveform_digis), "zsHPGe");} + event.put(std::move(ph_HPGe_digis), "phHPGe"); + + //LaBr + if (_saveRawWithHeaderWaveform_LaBr){ event.put(std::move(raw_LaBr_header_waveform_digis), "rawWithHeaderLaBr"); } + if (_saveRawWaveform_LaBr){event.put(std::move(raw_LaBr_waveform_digis), "rawLaBr");} + if (_saveZSWaveform_LaBr){event.put(std::move(zs_LaBr_waveform_digis), "zsLaBr");} + event.put(std::move(ph_LaBr_digis), "phLaBr"); + +} // produce() + +// ====================================================================== + + +void STMDigisFromFragments::endJob() { + if (_verbosityLevel >= 1){ + //Tells us what happened at the very end + std::cout << "\n========== STM JOB SUMMARY - (Unpacking Module) ==========\n"; + + std::cout << "Total events : " << _totalEvents << "\n"; + std::cout << "Total events w/HPGe&LaBr : " << _totalEventsWithHPGeLaBr << "\n"; + std::cout << "Total events w/only HPGe : " << _totalEventsWithOnlyHPGe << "\n"; + std::cout << "Total events w/only LaBr : " << _totalEventsWithOnlyLaBr << "\n"; + std::cout << "Total events w/Neither : " << _totalEventsWithNone << "\n"; + std::cout << "Total frags : " << _totalFragments << "\n"; + std::cout << "Total Container frags : " << _totalContainers << "\n"; + std::cout << "Total HPGe Containers : " << _totalContainersHPGe << "\n"; + std::cout << "Total LaBr Containers : " << _totalContainersLaBr << "\n"; + std::cout << "Total Inner frags : " << _totalInner << "\n"; + std::cout << "Unread frags : " << _unreadInnerFrags << "\n"; + + std::cout << "\n--- Data types read ---\n"; + std::cout << "RAW seen : " << _totalRaw << "\n"; + std::cout << "ZS seen : " << _totalZS << "\n"; + std::cout << "PH seen : " << _totalPH << "\n"; + std::cout << "\n"; + std::cout << "RAW LaBr : " << _totalRawLaBr << "\n"; + std::cout << "ZS LaBr : " << _totalZSLaBr << "\n"; + std::cout << "PH LaBr : " << _totalPHLaBr << "\n"; + std::cout << "\n"; + std::cout << "RAW HPGe : " << _totalRawHPGe << "\n"; + std::cout << "ZS HPGe : " << _totalZSHPGe << "\n"; + std::cout << "PH HPGe : " << _totalPHHPGe << "\n"; + + std::cout << "\n--- Data types filtered ---\n"; + std::cout << "Good RAW HPGe frags : " << _totalGoodRawHPGe << "\n"; + std::cout << "Good ZS HPGe frags : " << _totalGoodZSHPGe << "\n"; + std::cout << "Good PH HPGe frags : " << _totalGoodPHHPGe << "\n"; + std::cout << "\n"; + std::cout << "Zero RAW HPGe frags : " << _totalZeroRawHPGe << "\n"; + std::cout << "Zero ZS HPGe frags : " << _totalZeroZSHPGe << "\n"; + std::cout << "Zero PH HPGe frags : " << _totalZeroPHHPGe << "\n"; + std::cout << "\n"; + std::cout << "Empty Raw HPGe frags : " << _totalEmptyRawHPGe << "\n"; + std::cout << "Empty ZS HPGe frags : " << _totalEmptyZSHPGe << "\n"; + std::cout << "Empty PH HPGe frags : " << _totalEmptyPHHPGe << "\n"; + + std::cout << "\n"; + + std::cout << "Good RAW LaBr frags : " << _totalGoodRawLaBr << "\n"; + std::cout << "Good ZS LaBr frags : " << _totalGoodZSLaBr << "\n"; + std::cout << "Good PH LaBr frags : " << _totalGoodPHLaBr << "\n"; + std::cout << "\n"; + std::cout << "Zero RAW LaBr frags : " << _totalZeroRawLaBr << "\n"; + std::cout << "Zero ZS LaBr frags : " << _totalZeroZSLaBr << "\n"; + std::cout << "Zero PH LaBr frags : " << _totalZeroPHLaBr << "\n"; + std::cout << "\n"; + std::cout << "Empty Raw LaBr frags : " << _totalEmptyRawLaBr << "\n"; + std::cout << "Empty ZS LaBr frags : " << _totalEmptyZSLaBr << "\n"; + std::cout << "Empty PH LaBr frags : " << _totalEmptyPHLaBr << "\n"; + + std::cout << "=================================\n"; + + } +} + +DEFINE_ART_MODULE(STMDigisFromFragments) diff --git a/DAQ/src/STMPrintFragments_module.cc b/DAQ/src/STMPrintFragments_module.cc new file mode 100644 index 0000000000..6b833569a3 --- /dev/null +++ b/DAQ/src/STMPrintFragments_module.cc @@ -0,0 +1,97 @@ +// ====================================================================== +// +// STMPrintFragments_plugin: Add CRV data products to the event +// +// ====================================================================== + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "fhiclcpp/ParameterSet.h" + +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "art/Framework/Principal/Handle.h" +#include "artdaq-core-mu2e/Overlays/STMFragment.hh" +#include +#include "artdaq-core/Data/ContainerFragment.hh" + +#include +#include +#include + +namespace art +{ + class STMPrintFragments; +} + +using art::STMPrintFragments; + +// ====================================================================== + +class art::STMPrintFragments : public EDAnalyzer +{ + public: + struct Config + { + fhicl::Atom stmTag {fhicl::Name("stmTag"), fhicl::Comment("stmTag for new file")}; + }; + + // --- C'tor/d'tor: + explicit STMPrintFragments(const art::EDAnalyzer::Table& config); + + // --- Production: + virtual void analyze(const Event&); + + private: + + art::InputTag _stmFragmentsTag; +}; // STMPrintFragments + +// ====================================================================== + +STMPrintFragments::STMPrintFragments(const art::EDAnalyzer::Table& config) : + art::EDAnalyzer{config} + ,_stmFragmentsTag(config().stmTag()) +{} + +// ---------------------------------------------------------------------- + +void STMPrintFragments::analyze(const Event& event) +{ + art::EventNumber_t eventNumber = event.event(); + + auto STMContainerFragments = event.getValidHandle(_stmFragmentsTag); //Changed from STMFragments -> STMContainerFragments + + std::cout << std::dec << "Analyzer: Run " << event.run() << ", subrun " << event.subRun() + << ", event " << eventNumber << " has " << std::endl; + std::cout << STMContainerFragments->size() << " STM fragments." << std::endl; + + int frag_counter = 0; + for (auto& frag : *STMContainerFragments) { + ++frag_counter; + + //New lines + artdaq::ContainerFragment contf(frag); // interpret the fragment as a ContainerFragemnt (Will look inside here) + std::cout<<"N Blocks in the container = " << contf.block_count() << std::endl; //Should be 3 for the 3 STM Fragments + + for (size_t ii = 0; ii< contf.block_count(); ++ii){ + auto inner = contf.at(ii); + const auto dataBegin = inner->dataBegin(); + const auto dataEnd = inner->dataEnd(); + auto frag_id = inner->fragmentID(); + const auto stmDataBegin = reinterpret_cast(dataBegin); + const auto stmDataEnd = reinterpret_cast(dataEnd); + std::cout << "Frag_ID = " << frag_id << std::endl; + std::cout << "Container block_count = "< - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// ROOT includes -#include -#include -#include -#include -#include -#include -#include -#include - -namespace art -{ - class STMWaveformDigisFromFragments; -} - -using art::STMWaveformDigisFromFragments; - -// ====================================================================== - -class art::STMWaveformDigisFromFragments : public EDProducer -{ - public: - struct Config - { - fhicl::Atom stmTag {fhicl::Name("stmTag"), fhicl::Comment("Input module")}; - fhicl::Atom width_guess{ fhicl::Name("width_guess"), fhicl::Comment("Estimate of pulse width for fit (default = 150e-9)")}; - fhicl::Atom rise_time_guess{ fhicl::Name("rise_time_guess"), fhicl::Comment("Estimate of pulse rise time for fit (default = 6e-9)")}; - fhicl::Atom threshold{ fhicl::Name("threshold"), fhicl::Comment("Threshold to define the peak (default = 20000)")}; - fhicl::Atom samp_freq{ fhicl::Name("samp_freq"), fhicl::Comment("Sampling frequency of ADC (default = 300e6)")}; - fhicl::Atom pulse_region{ fhicl::Name("pulse_region"), fhicl::Comment("Number of ADC values to fit around pulse (default 51)")}; - fhicl::Atom pulses_per_event{fhicl::Name("pulses_per_event"), fhicl::Comment("Number of pulses per event to find Delta t (default = 2)")}; - }; - - // --- C'tor/d'tor: - explicit STMWaveformDigisFromFragments(const art::EDProducer::Table& config); - - // --- Production: - virtual void produce(Event&); - - private: - - art::InputTag _stmFragmentsTag; - - // Fhicl params - double _width_guess; - double _rise_time_guess; - double _threshold; - double _samp_freq; - int _pulse_region; - int _pulses_per_event; - - // Averaging params - // Pulse time average num - double avg_num = 0; - // Pulse time average den - double avg_den = 0; - - // Params to store end of last event - int lastEventLength = 0; - uint64_t lastEWT = 0; - std::vector prevEvent; - - int eventCounter = 0; - - static double rising_edge(double*t, double*p); - static double osc_rising_edge(double*t, double*p); - double fit_rising_edge(std::vector x, std::vector y, double pulse_time); - std::vector linspace(double start_in, double end_in, int num_in); - -}; // STMWaveformDigisFromFragments - -// ====================================================================== - -// Rising edge function -double STMWaveformDigisFromFragments::rising_edge(double* t, double* p){ - - double offset = p[0]; - double amplitude = p[1]; - double rise_time = p[2]; - double t0 = p[3]; - return offset + amplitude * (1 - exp(-(t[0] - t0) / rise_time)) * (t[0] >= t0); - -} - -// Fitting function -double STMWaveformDigisFromFragments::fit_rising_edge(std::vector x, std::vector y, double pulse_time) { - - // Function guess values - double offset_guess = *std::min_element(y.begin(),y.end()); // offset - double amp_guess = *std::max_element(y.begin(),y.end()); // amplitude - double center_guess = pulse_time; // pulse centre - double t0_guess = center_guess - _width_guess / 2; // rising edge t0 - - // Find the middle of the pulse (same as middle of data) - int middle = (x.size() / 2); - // x values to middle - std::vector fit_x(x.begin(), x.begin()+middle); - // y values to middle - std::vector fit_y(y.begin(), y.begin()+middle); - // Get start point - double start = fit_x.front(); - // Get end point - double end = fit_x.back(); - - int vsize = fit_x.size(); - - // Fit all pulses - TF1* fn = new TF1("fn", rising_edge, start, end, 4); - // Initial guesses for the parameters - fn->SetParameters(offset_guess, amp_guess, _rise_time_guess, t0_guess); - fn->SetParNames("Offset","A","#tau","t_{0}"); - fn->SetLineWidth(3); - - // Draw and fit graph - TGraph *g = new TGraph(vsize, &fit_x[0], &fit_y[0]); - g->Draw("A*"); - g->Fit("fn","Q"); - double rise_time = fn->GetParameter(2); - double t0 = fn->GetParameter(3); - // Calculate rising edge - double rising_edge_fit = t0 + rise_time; - // Return rising edge - return rising_edge_fit; -} - -std::vector STMWaveformDigisFromFragments::linspace(double start_in, double end_in, int num_in) -{ - - std::vector linspaced; - - double start = start_in; - double end = end_in; - double num = num_in; - - if (num == 0) { return linspaced; } - if (num == 1) - { - linspaced.push_back(start); - return linspaced; - } - - double delta = (end - start) / (num - 1); - - for(int i=0; i < num-1; ++i) - { - linspaced.push_back(start + delta * i); - } - linspaced.push_back(end); // I want to ensure that start and end - // are exactly the same as the input - return linspaced; -} - -STMWaveformDigisFromFragments::STMWaveformDigisFromFragments(const art::EDProducer::Table& config) : - art::EDProducer{config} - ,_stmFragmentsTag(config().stmTag()) - ,_width_guess(config().width_guess()) - ,_rise_time_guess(config().rise_time_guess()) - ,_threshold(config().threshold()) - ,_samp_freq(config().samp_freq()) - ,_pulse_region(config().pulse_region()) - ,_pulses_per_event(config().pulses_per_event()) - -{ - // Set the size of the vector - prevEvent.resize(_pulse_region); - produces(); -} - -// ---------------------------------------------------------------------- - - -void STMWaveformDigisFromFragments::produce(Event& event) -{ - eventCounter++; - std::unique_ptr stm_waveform_digis(new mu2e::STMWaveformDigiCollection); - - //auto STMFragments = event.getValidHandle("daq:STM"); - art::Handle STMFragmentsH; - if(!event.getByLabel(_stmFragmentsTag, STMFragmentsH)) { - event.put(std::move(stm_waveform_digis)); - return; - } - const auto STMFragments = STMFragmentsH.product(); - - int16_t detID = 0; - for (const auto& frag : *STMFragments) { - const auto& stm_frag = static_cast(frag); - - uint64_t EWT = *(stm_frag.EvNum()); - int16_t event_len = *(stm_frag.EvLen()); - detID = *(stm_frag.detID()) & 0xFF; - //std::cout << "DetID: " << *(stm_frag.detID()) << " | " << detID << "\n"; - //std::cout << "In event: " << event_number << " | Length of event: " << event_len << "\n"; - - // Get 40MHz DTC clock time - uint16_t DTC0 = static_cast((*(stm_frag.GetTHdr() + 28) >> 8) & 0xF) & 0xFFFF; - uint16_t DTC1 = static_cast(*(stm_frag.GetTHdr() + 29)) & 0xFFFF; - uint16_t DTC2 = static_cast(*(stm_frag.GetTHdr() + 30)) & 0xFFFF; - uint16_t DTC3 = static_cast(*(stm_frag.GetTHdr() + 31)) & 0xFFFF; - uint64_t DTC = static_cast(DTC3) << 40 | static_cast(DTC2) << 24 | static_cast(DTC1) << 8 | static_cast(DTC0); - //std::cout << "DTCs: " << std::hex << DTC0 << " " << DTC1 << " " << DTC2 << " " << DTC3 << " --> " << DTC << std::dec << std::endl; - //std::cout << "DTC time = " << DTC << " = " << (DTC/40e6) << "\n"; - - // Get 75MHz ADC clock time - int16_t adc[4] = {*(stm_frag.GetTHdr() + 4),*(stm_frag.GetTHdr() + 5),*(stm_frag.GetTHdr() + 6),*(stm_frag.GetTHdr() + 7)}; - uint64_t ADC; - memcpy(&ADC, adc, sizeof(ADC)); - //std::cout << "ADC time = " << ADC << " = " << (ADC/75e6) << "\n"; - - std::vector adcs; - adcs.reserve(event_len); - for (int16_t i_adc_sample = 0; i_adc_sample < event_len; ++i_adc_sample) { - adcs.emplace_back(*(stm_frag.DataBegin()+i_adc_sample)); - } - - // Number of pulses - int pulse_num = 0; - - // Initialise vectors for pulse locations and fit vals - std::vector pulse_locs; - std::vector fit_vals; - - // Get and fit the pulses - for(size_t i=0; i < adcs.size(); i++){ - // If the data is above the pulse threshold - if (adcs[i] > _threshold) { - // If the pulse started in the previous event - if (avg_den > 0 && i == 0) { - // Reset the average number... - avg_num = 0; - // Sum negative value from the previous event - for (int k = -avg_den; k < 0; ++k) { - avg_num += k; - } - } // End if pulse in previous event - // Average pulse time numerator - avg_num += i; - // Average pulse time denonimator - avg_den += 1; - } else { - // Else if not in a pulse region - if (avg_num != 0) { - // Calculate the pulse time average - int avg_time = static_cast(avg_num / avg_den); - // Store current pulse location - double pulse_loc = avg_time; - // Minimum pulse fit region - int pulse_min = static_cast(pulse_loc - (_pulse_region - 1) / 2); - // Maximum pulse fit region - int pulse_max = static_cast(pulse_loc + (_pulse_region - 1) / 2); - // Get diff - int pulse_diff = (pulse_max+1) - pulse_min; - // x values to fit - double p_min = pulse_min / _samp_freq; - double p_max = pulse_max / _samp_freq; - std::vector x(_pulse_region); - // Find linspace in pulse region - x = linspace(p_min, p_max, _pulse_region); - // y values to fit - std::vector y; - for(int k=0; k < pulse_diff; k++){ - y.push_back(adcs[pulse_min + k]); - } - // If the minimum pulse fit region is in the previous event - if (pulse_min < 0) { - // get start value in last event to copy - int del = pulse_min + _pulse_region; - // Store the y data that is in the previous event - std::vector prev_event(prevEvent.begin() + del, prevEvent.end()); - // Store the y data that is in this event - std::vector this_event(adcs.begin(), adcs.begin() + pulse_max + 1); - // New y values to fit - y.insert(y.begin(), prev_event.begin(), prev_event.end()); - y.insert(y.end(), this_event.begin(), this_event.end()); - } - // Fit rising edge - double rising_edge_fit = fit_rising_edge(x, y, pulse_loc / _samp_freq); - if(rising_edge_fit > 0){ - // rising edge is not in last event - // Save pulse time and fit vals in this event - pulse_locs.push_back(pulse_loc); - fit_vals.push_back(rising_edge_fit); - } - // Save the length of the event - lastEventLength = adcs.size(); - // Save last pulse region of event - int copyloc = adcs.size() - _pulse_region; - prevEvent.insert(prevEvent.begin(), adcs.begin() + copyloc, adcs.end()); - // Increment number of pulses - pulse_num++; - } - // Reset average variables - avg_num = 0; - avg_den = 0; - } - } - // Make sure variables are reset if events are not consecutive - if(EWT != lastEWT+1){ - avg_num = 0; - avg_den = 0; - } - - double peak_fitTime1 = 0; - double peak_fitTime2 = 0; - double peak_sep = 0; - // Find separation of first two consecutive pulses - if(pulse_locs.size()>2){ - peak_fitTime1 = fit_vals[0]*1e9; - peak_fitTime2 = fit_vals[1]*1e9; - peak_sep = (peak_fitTime2 - peak_fitTime1); - } - std::cout << "pulse 0 = " << peak_fitTime1 << " | pulse 1 = " << peak_fitTime2 << " | pulse sep = " << peak_sep << "\n"; - - // Store EWT - lastEWT = EWT; - - // Create the STMWaveformDigi and put it in the event - uint32_t trig_time_offset = 0; - mu2e::STMWaveformDigi stm_waveform(detID, EWT, DTC, ADC, trig_time_offset, peak_fitTime1, peak_fitTime2, peak_sep, adcs); - stm_waveform_digis->push_back(stm_waveform); - } - - event.put(std::move(stm_waveform_digis)); - -} // produce() - -// ====================================================================== - -DEFINE_ART_MODULE(STMWaveformDigisFromFragments) diff --git a/DAQ/test/inspectSTMFile.fcl b/DAQ/test/inspectSTMFile.fcl index 81d9e9737e..57f2fa8585 100644 --- a/DAQ/test/inspectSTMFile.fcl +++ b/DAQ/test/inspectSTMFile.fcl @@ -9,8 +9,6 @@ process_name : STMDigiFromFragment source : { module_type : RootInput fileNames : @nil - - # firstEvent : 500 maxEvents : -1 } @@ -18,39 +16,60 @@ source : { physics : { producers : { - makeSTMDigisHPGe : { - module_type : STMWaveformDigisFromFragments - stmTag : "daq:STM" - } + makeSTMDigis : { + module_type : STMDigisFromFragments + stmTag : "daq:ContainerSTM" + verbosityLevel : 1 + saveRawWithHeaderWaveform_HPGe: false + saveRawWaveform_HPGe: false + saveZSWaveform_HPGe: true + + saveRawWithHeaderWaveform_LaBr: false + saveRawWaveform_LaBr: false + saveZSWaveform_LaBr: true + } + + makeSTMBinary : { + module_type : STMBinaryDigisFromFragments + stmTag : "daq:ContainerSTM" + verbosityLevel : 1 + rawFile : "raw.bin" + zsFile : "zs.bin" + phFile : "ph.bin" + rawHeaderFile : "rawWithHeader.bin" + eventFile : "event.bin" + } } - analyzers : { - frag: - { - module_type : FragmentWatcher - mode_bitmask : 7 - } - #stm : - #{ - # module_type : STMPrintFragments - #} + analyzers : { + + stmPrint : { + module_type : STMPrintFragments + stmTag : "daq:ContainerSTM" + } } - t1 : [ makeSTMDigisHPGe ] - # e1 : [ stm, stmOutput ] - e1 : [ frag, stmOutput ] + t1 : [ makeSTMDigis ] + e1 : [ stmOutput ] + #e1 : [ ] + #e1 : [ frag, stmOutput ] trigger_paths : [t1] end_paths : [e1] } +# print summaries +services.scheduler.wantSummary: true +services.TimeTracker.printSummary: true + outputs : { stmOutput : { module_type : RootOutput fileName : "dig.stm.art" - outputCommands : [ "drop *_*_*_*", "keep mu2e::STMWaveformDigis_*_*_*" ] + outputCommands : [ "keep *_*_*_*", "drop artdaq::*_*_*_*" ] +# outputCommands : [ "drop *_*_*_*", "keep mu2e::STMWaveformDigis_*_*_*" ] } } diff --git a/RecoDataProducts/inc/STMMWDDigi.hh b/RecoDataProducts/inc/STMPHDigi.hh similarity index 58% rename from RecoDataProducts/inc/STMMWDDigi.hh rename to RecoDataProducts/inc/STMPHDigi.hh index 02cdce60df..d6d5fc9c1c 100644 --- a/RecoDataProducts/inc/STMMWDDigi.hh +++ b/RecoDataProducts/inc/STMPHDigi.hh @@ -1,5 +1,5 @@ -#ifndef RecoDataProducts_STMMWDDigi_hh -#define RecoDataProducts_STMMWDDigi_hh +#ifndef RecoDataProducts_STMPHDigi_hh +#define RecoDataProducts_STMPHDigi_hh // // Data product that represents an uncalibrated hit in the STM // @@ -12,11 +12,11 @@ namespace mu2e { - class STMMWDDigi { + class STMPHDigi { public: - STMMWDDigi() : _time(0), _energy(0){}; + STMPHDigi() : _time(0), _energy(0){}; - STMMWDDigi(uint32_t time, int16_t energy) : _time(time), _energy(energy) {}; + STMPHDigi(uint32_t time, int16_t energy) : _time(time), _energy(energy) {}; uint32_t time() const { return _time; } @@ -27,14 +27,14 @@ namespace mu2e { int16_t _energy; // uncalibrated energy [ADC units] }; - typedef std::vector STMMWDDigiCollection; + typedef std::vector STMPHDigiCollection; - bool lessByTime(const STMMWDDigi& a, const STMMWDDigi& b) { + bool lessByTime(const STMPHDigi& a, const STMPHDigi& b) { if (a.time() < b.time()) { return true; } else { return false; } } - bool lessByEnergy(const STMMWDDigi& a, const STMMWDDigi& b) { + bool lessByEnergy(const STMPHDigi& a, const STMPHDigi& b) { if (a.energy() < b.energy()) { return true; } else { return false; } } diff --git a/RecoDataProducts/inc/STMWaveformDigi.hh b/RecoDataProducts/inc/STMWaveformDigi.hh index 276b923df4..3c5f3fa003 100644 --- a/RecoDataProducts/inc/STMWaveformDigi.hh +++ b/RecoDataProducts/inc/STMWaveformDigi.hh @@ -14,40 +14,33 @@ #include "Offline/DataProducts/inc/STMChannel.hh" namespace mu2e { - class STMWaveformDigi { + public: - // Initialise all variables - STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(std::vector()) {}; - // Constructor for timing plus trig offset - STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _adcs(adcs) {}; - // Constructor for peak fitting - STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, double peak_fitTime1, double peak_fitTime2, double peak_sep, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(peak_fitTime1), _peak_fitTime2(peak_fitTime2), _peak_sep(peak_sep), _adcs(adcs) {}; - // Basic constructor - STMWaveformDigi(uint32_t trigTimeOffset, std::vector &adcs) : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(adcs) {}; + // Initialise all variables + STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _adcs(std::vector()) {}; + // Constructor for timing plus trig offset + STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), \ +_trigTimeOffset(trigTimeOffset), _adcs(adcs) {}; + // Basic constructor + STMWaveformDigi(uint32_t trigTimeOffset, std::vector &adcs) : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(trigTimeOffset), _adcs(adcs) {}; int16_t DetID () const { return _DetID; } uint64_t EWT () const { return _EWT; } uint64_t DTCtime() const { return _DTCtime; } uint64_t ADCtime() const { return _ADCtime; } uint32_t trigTimeOffset() const { return _trigTimeOffset; } - double peak_fitTime1 () const { return _peak_fitTime1; } - double peak_fitTime2 () const { return _peak_fitTime2; } - double peak_sep () const { return _peak_sep; } const std::vector& adcs () const { return _adcs; } - + void set_data ( size_t n_data, int16_t const* data ) { _adcs.resize(n_data); std::copy(data, data+n_data, _adcs.begin()); } + private: int16_t _DetID; uint64_t _EWT; uint64_t _DTCtime; uint64_t _ADCtime; - uint32_t _trigTimeOffset; // time offset from EWT? to first ADC value [ct] - double _peak_fitTime1; // fit time of first rising edge (ns) - double _peak_fitTime2; // fit time of second rising edge (ns) - double _peak_sep; // separation time (ns) - std::vector _adcs; // vector of ADC values for the waveform + uint32_t _trigTimeOffset; // time offset from EWT? to first ADC value [ct] + std::vector _adcs; // vector of ADC values for the waveform }; - typedef std::vector STMWaveformDigiCollection; } #endif diff --git a/RecoDataProducts/src/classes.h b/RecoDataProducts/src/classes.h index b95d8f6cd4..f191428da0 100644 --- a/RecoDataProducts/src/classes.h +++ b/RecoDataProducts/src/classes.h @@ -95,7 +95,7 @@ // STM #include "Offline/RecoDataProducts/inc/STMWaveformDigi.hh" -#include "Offline/RecoDataProducts/inc/STMMWDDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" #include "Offline/RecoDataProducts/inc/STMHit.hh" // MTP diff --git a/RecoDataProducts/src/classes_def.xml b/RecoDataProducts/src/classes_def.xml index 7e239ba0be..81380ec4f5 100644 --- a/RecoDataProducts/src/classes_def.xml +++ b/RecoDataProducts/src/classes_def.xml @@ -445,9 +445,9 @@ - - - + + + diff --git a/STMMC/CMakeLists.txt b/STMMC/CMakeLists.txt index 7e719143d7..ccec9dee51 100644 --- a/STMMC/CMakeLists.txt +++ b/STMMC/CMakeLists.txt @@ -43,8 +43,8 @@ cet_build_plugin(HPGeWaveformsFromStepPointMCs art::module Offline::STMConditions Offline::GlobalConstantsService ) -cet_build_plugin(MWDTree art::module - REG_SOURCE src/MWDTree_module.cc +cet_build_plugin(PHTree art::module + REG_SOURCE src/PHTree_module.cc LIBRARIES REG art_root_io::TFileService_service Offline::RecoDataProducts diff --git a/STMMC/fcl/HPGeReco.fcl b/STMMC/fcl/HPGeReco.fcl index b6dc9d5362..4627443482 100644 --- a/STMMC/fcl/HPGeReco.fcl +++ b/STMMC/fcl/HPGeReco.fcl @@ -49,9 +49,9 @@ physics : { makeTTreeGradients: false makeTTreeWaveforms: false } - MWDHPGe : { + PHHPGe : { module_type : STMMovingWindowDeconvolution - stmWaveformDigisTag : @local::STMMCAnalysis.MWD.HPGe.stmWaveformDigisTag.concatenated + stmWaveformDigisTag : @local::STMMCAnalysis.PH.HPGe.stmWaveformDigisTag.concatenated tau : @local::STMMCAnalysis.MWD.HPGe.tau M : @local::STMMCAnalysis.MWD.HPGe.M L : @local::STMMCAnalysis.MWD.HPGe.L @@ -59,7 +59,7 @@ physics : { thresholdgrad : @local::STMMCAnalysis.MWD.HPGe.thresholdgrad defaultBaselineMean : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineMean.suppressed defaultBaselineSD : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineSD.suppressed - makeTTreeMWD: false + makeTTreePH: false makeTTreeEnergies: false TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin verbosityLevel : 0 @@ -72,7 +72,7 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } - digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] + digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, PHHPGe] trigger_paths : [digitization_path] output_path : [compressedOutput] end_paths : [output_path] @@ -85,7 +85,7 @@ outputs : { SelectEvents: ["digitization_path"] outputCommands: [ "drop *_*_*_*", - "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" + "keep mu2e::STMPHDigis_PHHPGe_*_HPGeReco" ] } } diff --git a/STMMC/fcl/MakeTree.fcl b/STMMC/fcl/MakeTree.fcl index 1bcf442025..bad498198d 100644 --- a/STMMC/fcl/MakeTree.fcl +++ b/STMMC/fcl/MakeTree.fcl @@ -43,9 +43,9 @@ physics: { SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag Detector: @local::SimplifyStage2Data.DetectorName.LaBr } - MWDSpectra : { - module_type : MWDTree - STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag + PHSpectra : { + module_type : PHTree + STMPHDigiTag : @local::STMMCAnalysis.PH.HPGe.STMPHDigiTag EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin } } diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 1b56b4f7d1..2dca0648fa 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -111,7 +111,7 @@ STMMCAnalysis : { naverage : @local::STM.HPGe.naverage } } - MWD : { + PH : { HPGe : { stmWaveformDigisTag : { ZS: "ZSHPGe:" @@ -130,7 +130,7 @@ STMMCAnalysis : { full: @local::STM.HPGe.defaultBaselineSD suppressed: 0 } - STMMWDDigiTag : "MWDHPGe:" + STMPHDigiTag : "PHHPGe:" } LaBr : { stmWaveformDigisTag : "ZSLaBr:" @@ -141,7 +141,7 @@ STMMCAnalysis : { thresholdgrad : @local::STM.HPGe.thresholdgrad defaultBaselineMean : @local::STM.HPGe.defaultBaselineMean defaultBaselineSD : @local::STM.HPGe.defaultBaselineSD - STMMWDDigiTag : "MWDLaBr:" + STMPHDigiTag : "PHLaBr:" } } } diff --git a/STMMC/src/MWDTree_module.cc b/STMMC/src/PHTree_module.cc similarity index 72% rename from STMMC/src/MWDTree_module.cc rename to STMMC/src/PHTree_module.cc index e0703eca35..727192fb86 100644 --- a/STMMC/src/MWDTree_module.cc +++ b/STMMC/src/PHTree_module.cc @@ -1,5 +1,5 @@ // Adapted from ReadVirtualDetector_module.cc -// For STMMWDDigis, generates a TTree with the time and energy +// For STMPHDigis, generates a TTree with the time and energy // Original author: Ivan Logashenko // Adapted by: Pawel Plesniak @@ -25,7 +25,7 @@ #include "messagefacility/MessageLogger/MessageLogger.h" // Offline includes -#include "Offline/RecoDataProducts/inc/STMMWDDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" // ROOT includes #include "art_root_io/TFileService.h" @@ -36,44 +36,44 @@ typedef cet::map_vector_key key_type; typedef unsigned long VolumeId_type; namespace mu2e { - class MWDTree : public art::EDAnalyzer { + class PHTree : public art::EDAnalyzer { public: using Name=fhicl::Name; using Comment=fhicl::Comment; struct Config { - fhicl::Atom STMMWDDigiTag{Name("STMMWDDigiTag"), Comment("Tag identifying the MWD Digis")}; + fhicl::Atom STMPHDigiTag{Name("STMPHDigiTag"), Comment("Tag identifying the PH Digis")}; fhicl::OptionalAtom EnergyCalib{ Name("EnergyCalib"), Comment("Controls whether to make the energy TTrees with units of energy or in ADC values. If 0, will leave as ADC value, otherwise will multiply by this calibration to generate the energy.")}; }; using Parameters = art::EDAnalyzer::Table; - explicit MWDTree(const Parameters& conf); + explicit PHTree(const Parameters& conf); void analyze(const art::Event& e); void endJob(); private: - art::ProductToken STMMWDDigiToken; + art::ProductToken STMPHDigiToken; int eventCounter = 0, digiCounter = 0; TTree* ttree; uint32_t time = 0; double E = 0, EnergyCalib = 0; }; - MWDTree::MWDTree(const Parameters& conf) : + PHTree::PHTree(const Parameters& conf) : art::EDAnalyzer(conf), - STMMWDDigiToken(consumes(conf().STMMWDDigiTag())) { + STMPHDigiToken(consumes(conf().STMPHDigiTag())) { EnergyCalib = conf().EnergyCalib() ? *(conf().EnergyCalib()) : 1.0; art::ServiceHandle tfs; - ttree = tfs->make("ttree", "MWD ttree"); + ttree = tfs->make("ttree", "PH ttree"); ttree->Branch("time", &time, "time/i"); // ns ttree->Branch("E", &E, "E/D"); // keV }; - void MWDTree::analyze(const art::Event& event) { + void PHTree::analyze(const art::Event& event) { // Get the data products from the event - auto const& MWDDigis = event.getProduct(STMMWDDigiToken); - if (MWDDigis.empty()) + auto const& PHDigis = event.getProduct(STMPHDigiToken); + if (PHDigis.empty()) return; eventCounter++; // Loop over all VD hits - for (const STMMWDDigi& digi : MWDDigis) { + for (const STMPHDigi& digi : PHDigis) { // Extract the parameters time = digi.time(); E = digi.energy() * EnergyCalib; @@ -87,13 +87,13 @@ namespace mu2e { return; }; - void MWDTree::endJob() { - mf::LogInfo log("MWD tree summary"); - log << "=========MWD tree summary =========\n"; + void PHTree::endJob() { + mf::LogInfo log("PH tree summary"); + log << "========= PH tree summary =========\n"; log << std::left << std::setw(25) << "\tProcessed events: " << eventCounter << "\n"; log << std::left << std::setw(25) << "\tProcessed digis: " << digiCounter << "\n"; log << "===================================\n"; }; }; // end namespace mu2e -DEFINE_ART_MODULE(mu2e::MWDTree) +DEFINE_ART_MODULE(mu2e::PHTree) diff --git a/STMReco/CMakeLists.txt b/STMReco/CMakeLists.txt index 6488dcaa77..be76a02e26 100644 --- a/STMReco/CMakeLists.txt +++ b/STMReco/CMakeLists.txt @@ -17,8 +17,8 @@ cet_build_plugin(PlotSTMEnergySpectrum art::module Offline::RecoDataProducts ) -cet_build_plugin(PlotSTMMWDSpectrum art::module - REG_SOURCE src/PlotSTMMWDSpectrum_module.cc +cet_build_plugin(PlotSTMPHSpectrum art::module + REG_SOURCE src/PlotSTMPHSpectrum_module.cc LIBRARIES REG art_root_io::TFileService_service Offline::GlobalConstantsService @@ -71,7 +71,7 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/makeSTMHits_testbeam.fcl ${CURR configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/mwd.fcl ${CURRENT_BINARY_DIR} fcl/mwd.fcl COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/mwd_testbeam.fcl ${CURRENT_BINARY_DIR} fcl/mwd_testbeam.fcl COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/plotSTMEnergySpectrum.fcl ${CURRENT_BINARY_DIR} fcl/plotSTMEnergySpectrum.fcl COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/plotSTMMWDSpectrum.fcl ${CURRENT_BINARY_DIR} fcl/plotSTMMWDSpectrum.fcl COPYONLY) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/plotSTMPHSpectrum.fcl ${CURRENT_BINARY_DIR} fcl/plotSTMPHSpectrum.fcl COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/plotSTMWaveformDigis.fcl ${CURRENT_BINARY_DIR} fcl/plotSTMWaveformDigis.fcl COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/prolog.fcl ${CURRENT_BINARY_DIR} fcl/prolog.fcl COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/prolog_testbeam.fcl ${CURRENT_BINARY_DIR} fcl/prolog_testbeam.fcl COPYONLY) diff --git a/STMReco/fcl/makeSTMHits.fcl b/STMReco/fcl/makeSTMHits.fcl index 6785ddf8b6..78504400b9 100644 --- a/STMReco/fcl/makeSTMHits.fcl +++ b/STMReco/fcl/makeSTMHits.fcl @@ -1,5 +1,5 @@ # -# Create STMHits from STMMWDDigis +# Create STMHits from STMPHDigis # #include "Offline/fcl/standardServices.fcl" @@ -14,10 +14,10 @@ services : @local::Services.Reco physics: { producers : { stmHitsHPGe : { module_type : MakeSTMHits - stmMWDDigisTag : "mwdHPGe" + stmPHDigisTag : "makeSTMDigis:phHPGe" } stmHitsLaBr : { module_type : MakeSTMHits - stmMWDDigisTag : "mwdLaBr" + stmPHDigisTag : "makeSTMDigis:phLaBr" } } filters : { @@ -27,7 +27,7 @@ physics: { # setup paths HPGePath : [ stmHitsHPGe ] LaBrPath : [ stmHitsLaBr ] - trigger_paths: [ HPGePath, LaBrPath ] + trigger_paths: [LaBrPath ] outPath : [ STMOutput ] end_paths: [outPath] } diff --git a/STMReco/fcl/mwd.fcl b/STMReco/fcl/mwd.fcl index c7021aee3c..ddef0c8f45 100644 --- a/STMReco/fcl/mwd.fcl +++ b/STMReco/fcl/mwd.fcl @@ -5,7 +5,7 @@ #include "Offline/fcl/standardServices.fcl" #include "Offline/STMReco/fcl/prolog.fcl" -process_name: STMMWD +process_name: STMPH source : { module_type : RootInput @@ -17,34 +17,40 @@ services : { physics: { producers : { - mwdHPGe : { + phHPGe : { module_type : STMMovingWindowDeconvolution - stmWaveformDigisTag : "zeroSuppressHPGe" - verbosityLevel : 1 + stmWaveformDigisTag : "makeSTMDigis:zsHPGe" + verbosityLevel : 6 tau : @local::STM.HPGe.tau M : @local::STM.HPGe.M L : @local::STM.HPGe.L nsigma_cut : @local::STM.HPGe.nsigma_cut thresholdgrad : @local::STM.HPGe.thresholdgrad + defaultBaselineMean: @local::STM.HPGe.defaultBaselineMean + defaultBaselineSD: @local::STM.HPGe.defaultBaselineSD + xAxis: "sample_number" } - mwdLaBr : { + phLaBr : { module_type : STMMovingWindowDeconvolution - stmWaveformDigisTag : "zeroSuppressLaBr" - verbosityLevel : 1 + stmWaveformDigisTag : "makeSTMDigis:LaBr" + verbosityLevel : 6 tau : @local::STM.LaBr.tau M : @local::STM.LaBr.M L : @local::STM.LaBr.L nsigma_cut : @local::STM.LaBr.nsigma_cut thresholdgrad : @local::STM.LaBr.thresholdgrad + defaultBaselineMean: @local::STM.HPGe.defaultBaselineMean + defaultBaselineSD: @local::STM.HPGe.defaultBaselineSD + xAxis: "sample_number" } } filters : { } analyzers : { } # setup paths - HPGePath : [ mwdHPGe ] - LaBrPath : [ mwdLaBr ] + HPGePath : [ phHPGe ] + LaBrPath : [ phLaBr ] trigger_paths: [ HPGePath, LaBrPath ] outPath : [ STMOutput ] end_paths: [outPath] @@ -56,6 +62,6 @@ outputs: { module_type: RootOutput outputCommands: [ "keep *_*_*_*" ] SelectEvents: [ ] - fileName : "dig.owner.MWDWaveformDigis.version.sequencer.art" + fileName : "dig.owner.PHWaveformDigis.version.sequencer.art" } } diff --git a/STMReco/fcl/mwd_testbeam.fcl b/STMReco/fcl/mwd_testbeam.fcl index 3369288552..e2ca4fa4b7 100644 --- a/STMReco/fcl/mwd_testbeam.fcl +++ b/STMReco/fcl/mwd_testbeam.fcl @@ -1,5 +1,5 @@ # -# Apply mpoving window deconvolution algorithm to zero-suppressed STMWaveformDigis +# Apply Moving Window Deconvolution algorithm to zero-suppressed STMWaveformDigis # #include "Offline/STMReco/fcl/prolog_testbeam.fcl" diff --git a/STMReco/fcl/plotSTMMWDSpectrum.fcl b/STMReco/fcl/plotSTMMWDSpectrum.fcl deleted file mode 100644 index 3ed41cc3e6..0000000000 --- a/STMReco/fcl/plotSTMMWDSpectrum.fcl +++ /dev/null @@ -1,35 +0,0 @@ -# -# Plot STMWaveformDigis -# - -#include "Offline/fcl/standardServices.fcl" - -process_name: PlotSTMMWDSpectrum - -source : { - module_type : RootInput -} - -services : @local::Services.Reco -physics: { - producers : { - } - filters : { - } - analyzers : { - plotHPGeMWDSpectrum : { - module_type : PlotSTMMWDSpectrum - stmMWDDigisTag : "mwdHPGe" - } - plotLaBrMWDSpectrum : { - module_type : PlotSTMMWDSpectrum - stmMWDDigisTag : "mwdLaBr" - } - } - # setup paths - trigger_paths: [ ] - anaPath : [ plotHPGeMWDSpectrum, plotLaBrMWDSpectrum ] - end_paths: [anaPath] -} - -services.TFileService.fileName : "stmMWDSpectrum.root" diff --git a/STMReco/fcl/plotSTMPHSpectrum.fcl b/STMReco/fcl/plotSTMPHSpectrum.fcl new file mode 100644 index 0000000000..b38b70cec3 --- /dev/null +++ b/STMReco/fcl/plotSTMPHSpectrum.fcl @@ -0,0 +1,35 @@ +# +# Plot STMWaveformDigis +# + +#include "Offline/fcl/standardServices.fcl" + +process_name: PlotSTMPHSpectrum + +source : { + module_type : RootInput +} + +services : @local::Services.Reco +physics: { + producers : { + } + filters : { + } + analyzers : { + plotHPGePHSpectrum : { + module_type : PlotSTMPHSpectrum + stmPHDigisTag : "phHPGe" + } + plotLaBrPHSpectrum : { + module_type : PlotSTMPHSpectrum + stmPHDigisTag : "makeSTMDigis:phHPGe" + } + } + # setup paths + trigger_paths: [ ] + anaPath : [ plotLaBrPHSpectrum ] + end_paths: [anaPath] +} + +services.TFileService.fileName : "stmPHSpectrum.root" diff --git a/STMReco/fcl/plotSTMWaveformDigis.fcl b/STMReco/fcl/plotSTMWaveformDigis.fcl index f51f63401f..e6b1bca596 100644 --- a/STMReco/fcl/plotSTMWaveformDigis.fcl +++ b/STMReco/fcl/plotSTMWaveformDigis.fcl @@ -21,40 +21,59 @@ physics: { producers : { } filters : { } analyzers : { - plotHPGeWaveformDigis : { + + plotRawWaveformDigisHPGe : { module_type : PlotSTMWaveformDigis - stmWaveformDigisTag : "FromSTMTestBeamData:HPGe" + stmWaveformDigisTag : "makeSTMDigis:rawHPGe" subtractPedestal : false verbosityLevel : 0 - xAxis : "event_time" + xAxis : "sample_number" } - plotLaBrWaveformDigis : { + + plotZSWaveformDigisHPGe : { module_type : PlotSTMWaveformDigis - stmWaveformDigisTag : "FromSTMTestBeamData:LaBr" + stmWaveformDigisTag : "makeSTMDigis:zsHPGe" subtractPedestal : false verbosityLevel : 0 - xAxis : "event_time" + xAxis : "sample_number" } - plotHPGeWaveformDigisZP : { + plotPHWaveformDigisHPGe : { module_type : PlotSTMWaveformDigis - stmWaveformDigisTag : "zeroSuppressHPGe" + stmWaveformDigisTag : "makeSTMDigis:phHPGe" subtractPedestal : false verbosityLevel : 0 - xAxis : "event_time" + xAxis : "sample_number" } - plotLaBrWaveformDigisZP : { + + plotRawWaveformDigisLaBr : { + module_type : PlotSTMWaveformDigis + stmWaveformDigisTag : "makeSTMDigis:rawLaBr" + subtractPedestal : false + verbosityLevel : 0 + xAxis : "sample_number" + } + + plotZSWaveformDigisLaBr : { module_type : PlotSTMWaveformDigis - stmWaveformDigisTag : "zeroSuppressLaBr" + stmWaveformDigisTag : "makeSTMDigis:zsLaBr" subtractPedestal : false verbosityLevel : 0 - xAxis : "event_time" + xAxis : "sample_number" } + + plotPHWaveformDigisLaBr : { + module_type : PlotSTMWaveformDigis + stmWaveformDigisTag : "makeSTMDigis:phLaBr" + subtractPedestal : false + verbosityLevel : 0 + xAxis : "sample_number" + } + } # setup paths trigger_paths: [ ] - anaPath : [ plotHPGeWaveformDigis, plotLaBrWaveformDigis, - plotHPGeWaveformDigisZP, plotLaBrWaveformDigisZP ] + anaPath : [ plotRawWaveformDigisHPGe, plotZSWaveformDigisHPGe, plotRawWaveformDigisLaBr, plotZSWaveformDigisLaBr ] end_paths: [anaPath] } diff --git a/STMReco/src/MakeSTMHits_module.cc b/STMReco/src/MakeSTMHits_module.cc index 3336b4de69..abde7f2e4b 100644 --- a/STMReco/src/MakeSTMHits_module.cc +++ b/STMReco/src/MakeSTMHits_module.cc @@ -1,5 +1,5 @@ // -// Create STMHits from STMMWDDigis +// Create STMHits from STMPHDigis // #include "art/Framework/Principal/Event.h" #include "art/Framework/Core/EDProducer.h" @@ -22,7 +22,7 @@ #include "TH1F.h" #include "TTree.h" -#include "Offline/RecoDataProducts/inc/STMMWDDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" #include "Offline/RecoDataProducts/inc/STMHit.hh" // C++ @@ -37,7 +37,7 @@ namespace mu2e { using Name=fhicl::Name; using Comment=fhicl::Comment; struct Config { - fhicl::Atom stmMWDDigisTag{ Name("stmMWDDigisTag"), Comment("InputTag for STMMWDDigiCollection")}; + fhicl::Atom stmPHDigisTag{ Name("stmPHDigisTag"), Comment("InputTag for STMPHDigiCollection")}; }; using Parameters = art::EDProducer::Table; explicit MakeSTMHits(const Parameters& conf); @@ -45,33 +45,34 @@ namespace mu2e { private: void produce(art::Event& e) override; - art::ProductToken _stmMWDDigisToken; + art::ProductToken _stmPHDigisToken; STMChannel _channel; ProditionsHandle _stmEnergyCalib_h; }; MakeSTMHits::MakeSTMHits(const Parameters& config ) : art::EDProducer{config} - ,_stmMWDDigisToken(consumes(config().stmMWDDigisTag())) - ,_channel(STMUtils::getChannel(config().stmMWDDigisTag())) + ,_stmPHDigisToken(consumes(config().stmPHDigisTag())) + ,_channel(STMUtils::getChannel(config().stmPHDigisTag())) ,_stmEnergyCalib_h() { produces(); } - + //Originally had _channel(STMChannel::LaBr) + void MakeSTMHits::produce(art::Event& event) { // create output unique_ptr outputSTMHits(new STMHitCollection); - auto mwdDigisHandle = event.getValidHandle(_stmMWDDigisToken); + auto phDigisHandle = event.getValidHandle(_stmPHDigisToken); STMEnergyCalib const& stmEnergyCalib = _stmEnergyCalib_h.get(event.id()); // get calibration const auto nsPerCt = stmEnergyCalib.nsPerCt(_channel); const auto& pars = stmEnergyCalib.calib(_channel); - for (const auto& mwd_digi : *mwdDigisHandle) { - auto uncalib_time = mwd_digi.time(); - auto uncalib_energy = mwd_digi.energy(); + for (const auto& ph_digi : *phDigisHandle) { + auto uncalib_time = ph_digi.time(); + auto uncalib_energy = ph_digi.energy(); float time = uncalib_time*nsPerCt; float energy = pars.p0 + pars.p1*uncalib_energy + pars.p2*uncalib_energy*uncalib_energy; diff --git a/STMReco/src/PlotSTMEnergySpectrum_module.cc b/STMReco/src/PlotSTMEnergySpectrum_module.cc index 1a0fa43424..64850d1e83 100644 --- a/STMReco/src/PlotSTMEnergySpectrum_module.cc +++ b/STMReco/src/PlotSTMEnergySpectrum_module.cc @@ -1,5 +1,5 @@ // -// Analyzer module to create a histogram of the STMMWDDigi uncalibrated energies +// Analyzer module to create a histogram of the STMPHDigi uncalibrated energies // #include "art/Framework/Principal/Event.h" #include "art/Framework/Core/EDAnalyzer.h" diff --git a/STMReco/src/PlotSTMMWDSpectrum_module.cc b/STMReco/src/PlotSTMMWDSpectrum_module.cc deleted file mode 100644 index e75051fc64..0000000000 --- a/STMReco/src/PlotSTMMWDSpectrum_module.cc +++ /dev/null @@ -1,72 +0,0 @@ -// -// Analyzer module to create a histogram of the STMMWDDigi uncalibrated energies -// -#include "art/Framework/Principal/Event.h" -#include "art/Framework/Core/EDAnalyzer.h" -#include "art/Framework/Principal/Handle.h" -#include "art/Framework/Core/ModuleMacros.h" -#include "cetlib_except/exception.h" -#include "fhiclcpp/types/Atom.h" -#include "canvas/Utilities/InputTag.h" -#include "messagefacility/MessageLogger/MessageLogger.h" - -#include "art_root_io/TFileService.h" -#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" -#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" - -#include "Offline/MCDataProducts/inc/StepPointMC.hh" -#include -// root -#include "TH1F.h" -#include "TF1.h" -#include "TTree.h" -#include "TSpectrum.h" -#include "TGraph.h" - -#include "Offline/RecoDataProducts/inc/STMMWDDigi.hh" - -using namespace std; -using CLHEP::Hep3Vector; -namespace mu2e { - - class PlotSTMMWDSpectrum : public art::EDAnalyzer { - public: - using Name=fhicl::Name; - using Comment=fhicl::Comment; - struct Config { - fhicl::Atom stmMWDDigisTag{ Name("stmMWDDigisTag"), Comment("InputTag for STMMWDDigiCollection")}; - }; - using Parameters = art::EDAnalyzer::Table; - explicit PlotSTMMWDSpectrum(const Parameters& conf); - - private: - void beginJob() override; - void analyze(const art::Event& e) override; - - TH1D* _mwdSpectrum; - art::ProductToken _stmMWDDigisToken; - }; - - PlotSTMMWDSpectrum::PlotSTMMWDSpectrum(const Parameters& config ) : - art::EDAnalyzer{config}, - _stmMWDDigisToken(consumes(config().stmMWDDigisTag())) - { } - - void PlotSTMMWDSpectrum::beginJob() { - art::ServiceHandle tfs; - // create histograms - _mwdSpectrum=tfs->make("mwdSpectrum", "MWD Spectrum", 1000, 0, 1e4); - } - - void PlotSTMMWDSpectrum::analyze(const art::Event& event) { - - auto mwdDigisHandle = event.getValidHandle(_stmMWDDigisToken); - - for (const auto& mwdDigi : *mwdDigisHandle) { - auto energy = mwdDigi.energy(); - _mwdSpectrum->Fill(energy); - } - } -} - -DEFINE_ART_MODULE(mu2e::PlotSTMMWDSpectrum) diff --git a/STMReco/src/PlotSTMPHSpectrum_module.cc b/STMReco/src/PlotSTMPHSpectrum_module.cc new file mode 100644 index 0000000000..4995dc4aef --- /dev/null +++ b/STMReco/src/PlotSTMPHSpectrum_module.cc @@ -0,0 +1,82 @@ +// +// Analyzer module to create a histogram of the STMPHDigi uncalibrated energies +// +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "cetlib_except/exception.h" +#include "fhiclcpp/types/Atom.h" +#include "canvas/Utilities/InputTag.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "art_root_io/TFileService.h" +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" + +#include "Offline/MCDataProducts/inc/StepPointMC.hh" +#include +// root +#include "TH2F.h" +#include "TH1F.h" +#include "TF1.h" +#include "TTree.h" +#include "TSpectrum.h" +#include "TGraph.h" + +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" + +using namespace std; +using CLHEP::Hep3Vector; +namespace mu2e { + + class PlotSTMPHSpectrum : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom stmPHDigisTag{ Name("stmPHDigisTag"), Comment("InputTag for STMPHDigiCollection")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit PlotSTMPHSpectrum(const Parameters& conf); + + private: + void beginJob() override; + void analyze(const art::Event& e) override; + + TH2F* _twoDhist; //Histograms of Energy vs binned event + int eventCount = 0; + + TH1D* _phSpectrum; + art::ProductToken _stmPHDigisToken; + }; + + PlotSTMPHSpectrum::PlotSTMPHSpectrum(const Parameters& config ) : + art::EDAnalyzer{config}, + _stmPHDigisToken(consumes(config().stmPHDigisTag())) + { } + + void PlotSTMPHSpectrum::beginJob() { + art::ServiceHandle tfs; + // create histograms + _phSpectrum=tfs->make("phSpectrum", "PH Spectrum", 1000, 0, 1e4); + _twoDhist=tfs->make("twoDhist","Pulse Height vs events;Event Bins; Pulse Height", + 1000,0,1000, // X-axis scale + 1000,0,1e5); // Y-axis scale + } + + void PlotSTMPHSpectrum::analyze(const art::Event& event) { + + auto phDigisHandle = event.getValidHandle(_stmPHDigisToken); + int binBlock = eventCount/100; + + for (const auto& phDigi : *phDigisHandle) { + auto energy = phDigi.energy(); + _phSpectrum->Fill(energy); + _twoDhist->Fill(binBlock, energy); + } + ++eventCount; + } +} + +DEFINE_ART_MODULE(mu2e::PlotSTMPHSpectrum) diff --git a/STMReco/src/PlotSTMWaveformDigis_module.cc b/STMReco/src/PlotSTMWaveformDigis_module.cc index ba314322bc..84c1556bb2 100644 --- a/STMReco/src/PlotSTMWaveformDigis_module.cc +++ b/STMReco/src/PlotSTMWaveformDigis_module.cc @@ -9,11 +9,14 @@ #include "fhiclcpp/types/Atom.h" #include "canvas/Utilities/InputTag.h" #include "messagefacility/MessageLogger/MessageLogger.h" - #include "art_root_io/TFileService.h" #include #include +#include +#include +#include + // root #include "TH1F.h" #include "TF1.h" @@ -34,14 +37,20 @@ namespace mu2e { struct Config { fhicl::Atom stmWaveformDigisTag{ Name("stmWaveformDigisTag"), Comment("InputTag for STMWaveformDigiCollection")}; fhicl::Atom subtractPedestal{ Name("subtractPedestal"), Comment("True/False whether to subtract the pedestal before plotting")}; - fhicl::Atom xAxis{ Name("xAxis"), Comment("Choice of x-axis unit: \"sample_number\", \"adcs_time\", or \"event_time\"") }; + fhicl::Atom xAxis{ Name("xAxis"), Comment("Choice of x-axis unit: \"sample_number\", \"adcs_time\", or \"event_time\"")} ; fhicl::Atom verbosityLevel{ Name("verbosityLevel"), Comment("Verbosity level")}; }; using Parameters = art::EDAnalyzer::Table; explicit PlotSTMWaveformDigis(const Parameters& conf); - + private: + void beginJob() override;//For _hist + void endJob() override; //For printing counter void analyze(const art::Event& e) override; + + TH1F* _hist; //Hist for WaveLength + int _zeroLengthCount = 0; + art::InputTag _stmWaveformDigisTag; art::ProductToken _stmWaveformDigisToken; bool _subtractPedestal; @@ -53,18 +62,43 @@ namespace mu2e { PlotSTMWaveformDigis::PlotSTMWaveformDigis(const Parameters& config ) : art::EDAnalyzer{config}, + _stmWaveformDigisTag{config().stmWaveformDigisTag()}, _stmWaveformDigisToken(consumes(config().stmWaveformDigisTag())), _subtractPedestal(config().subtractPedestal()), _xAxis(config().xAxis()), _verbosityLevel(config().verbosityLevel()), + //_channel(STMChannel::findByName("HPGe")) _channel(STMUtils::getChannel(config().stmWaveformDigisTag())) + { } + void PlotSTMWaveformDigis::beginJob(){ + + art::ServiceHandle tfs; + std::string X = std::string(_stmWaveformDigisTag.instance()); //Gets instance name fromf fcl + std::transform(X.begin(),X.end(),X.begin(), toupper); //Raises uppercase of DigiTag + std::string hWaveLength_title = "Waveform Lengths for " + X + " Pulses"; //Builds title + _hist = tfs->make("hWaveLength", hWaveLength_title.c_str() ,1000,0,1000); //makes the histogram hWaveLength + + } + + void PlotSTMWaveformDigis::endJob(){ + if (_verbosityLevel > 1){ + std::cout << " Zero length Waveforms Count = " << _zeroLengthCount << std::endl; + } + } + void PlotSTMWaveformDigis::analyze(const art::Event& event) { + //Boolean for whether we get a match for zsHPGe or LaBr + const std::string instance = std::string(_stmWaveformDigisTag.instance()); + const bool plotZSOffsetWaveforms = (instance == "zsHPGe" || instance == "zsLaBr");// Don't need offset waveforms for raw waveformdigis + + if (!plotZSOffsetWaveforms){ + if (_verbosityLevel >1){ std::cout << "Instance : " << instance << " , not ZS ==> No Offset Waveform will be created" << std::endl;} } + art::ServiceHandle tfs; auto waveformsHandle = event.getValidHandle(_stmWaveformDigisToken); - std::stringstream histname, histtitle; int count = 0; STMEnergyCalib const& stmEnergyCalib = _stmEnergyCalib_h.get(event.id()); // get prodition @@ -74,29 +108,61 @@ namespace mu2e { } const auto nsPerCt = stmEnergyCalib.nsPerCt(_channel); + if (_verbosityLevel > 1){ std::cout<<"size = "<size()<make(histname.str().c_str(), histtitle.str().c_str(), binning.nbins(),binning.low(),binning.high()); - - for (size_t i_adc = 0; i_adc < adcs.adcs().size(); ++i_adc) { - const auto adc = adcs.adcs().at(i_adc); - - auto content = adc; - if (_subtractPedestal) { - content -= pedestal; - } - - hWaveform->SetBinContent(i_adc+1, content); // bins start numbering at 1 - } - ++count; - } - } + if (waveform.adcs().size() == 0){ + ++_zeroLengthCount; + } else { + //None empty waveforms go in here + _hist->Fill(waveform.adcs().size()); //_hist was created outside so there should be no problem here + + Binning binning = STMUtils::getBinning(waveform, _xAxis, nsPerCt); + TH1F* hWaveform = tfs->make(histname.str().c_str(), histtitle.str().c_str(), binning.nbins(),binning.low(),binning.high()); + TH1F* hWaveformOffset = nullptr; // Standby + + hWaveform->GetYaxis()->SetTitle("ADCs"); + hWaveform->GetXaxis()->SetTitle("Sample Number"); + + if (plotZSOffsetWaveforms){ + //For offset waveforms -> Better organize this area + const auto zs_offset = waveform.trigTimeOffset(); //Grabs stored offset + std::stringstream histname2;//New histogram to include the offset waveform from trigTimeOffset + std::stringstream histtitle2; + histname2 << histname.str() << "_offset" << zs_offset; + histtitle2 << histtitle2.str() << instance << " , offset : " << zs_offset; + + hWaveformOffset = tfs->make( histname2.str().c_str(), histtitle.str().c_str(), binning.nbins(), binning.low()+zs_offset, binning.high()+zs_offset );//Shifting bins using offset + hWaveformOffset->GetYaxis()->SetTitle("ADCs"); + hWaveformOffset->GetXaxis()->SetTitle("Sample Number (Includes + trigtimeOffset"); + // int n_bins = hWaveformOffset->GetNbinsX(); //Grabs already contained nbins from waveform + //double xmin = hWaveformOffset->GetXaxis()->GetXmin();// gets xmin from waveform + //double xmax = hWaveformOffset->GetXaxis()->GetXmax();//gets xman from waveform + //hWaveformOffset->SetBins(n_bins, xmin + zs_offset, xmax + zs_offset);// shifts the xmin and xmax by offset, keeps numbers of bins + }//PlotZSOffset + + for (size_t i_adc = 0; i_adc < waveform.adcs().size(); ++i_adc){ + + const auto adc = waveform.adcs().at(i_adc); + auto content = adc; + if (_subtractPedestal) { + content -= pedestal; + } + + hWaveform->SetBinContent(i_adc+1,content); + if ( plotZSOffsetWaveforms) { hWaveformOffset->SetBinContent(i_adc+1,content);} + + } + }//else + ++count; + }//for + }//analyzer } DEFINE_ART_MODULE(mu2e::PlotSTMWaveformDigis) diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 6da8c05fc4..73969557ae 100644 --- a/STMReco/src/STMMovingWindowDeconvolution_module.cc +++ b/STMReco/src/STMMovingWindowDeconvolution_module.cc @@ -35,7 +35,7 @@ #include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" #include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" #include "Offline/RecoDataProducts/inc/STMWaveformDigi.hh" -#include "Offline/RecoDataProducts/inc/STMMWDDigi.hh" +#include "Offline/RecoDataProducts/inc/STMPHDigi.hh" #include "Offline/Mu2eUtilities/inc/STMUtils.hh" #include "Offline/ProditionsService/inc/ProditionsHandle.hh" #include "Offline/STMConditions/inc/STMEnergyCalib.hh" @@ -61,11 +61,12 @@ namespace mu2e { fhicl::Atom thresholdgrad{Name("thresholdgrad"), Comment("Threshold on gradient to cut out peaks when calculating baseline")}; fhicl::Atom defaultBaselineMean{Name("defaultBaselineMean"), Comment("Default mean to use for baseline when ZS removes all baseline data, in ADC values")}; fhicl::Atom defaultBaselineSD{Name("defaultBaselineSD"), Comment("Default standard deviation to use for baseline when ZS removes all baseline data, in ADC values")}; - fhicl::OptionalAtom makeTTreeMWD{ Name("makeTTreeMWD"), Comment("Controls whether to make the TTree with branches ADC, deconvoluted, differentiated, averaged")}; + fhicl::OptionalAtom makeTTreePH{ Name("makeTTreePH"), Comment("Controls whether to make the TTree with branches ADC, deconvoluted, differentiated, averaged")}; fhicl::OptionalAtom makeTTreeEnergies{ Name("makeTTreeEnergies"), Comment("Controls whether to make the TTree with branches time, E")}; fhicl::OptionalAtom TTreeEnergyCalib{ Name("TTreeEnergyCalib"), Comment("Controls whether to make the energy TTrees with units of energy or in ADC values. If 0, will leave as ADC value, otherwise will multiply by this calibration to generate the energy.")}; fhicl::OptionalAtom verbosityLevel{Name("verbosityLevel"), Comment("Verbosity level")}; - fhicl::OptionalAtom xAxis{ Name("xAxis"), Comment("Choice of x-axis unit for histograms if verbosity level >= 5: \"sample_number\", \"waveform_time\", or \"event_time\"") }; + fhicl::OptionalAtom xAxis{ Name("xAxis"), Comment("Choice of x-axis unit for histograms if verbosity level >= 5: \"sample_number\", \"waveform_time\", or \"event_time\"") + }; }; using Parameters = art::EDProducer::Table; @@ -94,7 +95,7 @@ namespace mu2e { double thresholdgrad = 0.0; // threshold on gradient double defaultBaselineMean = 0.0; double defaultBaselineSD = 0.0; - bool makeTTreeMWD = false; // controls whether to make MWD process TTree + bool makeTTreePH = false; // controls whether to make PH process TTree bool makeTTreeEnergies = false; // controls whether to make results TTree int verbosityLevel = 0; // std cout verbosity level std::string _xAxis = ""; // optional parameter for x-axis unit if plotting histograms @@ -102,7 +103,7 @@ namespace mu2e { // Proditions service ProditionsHandle _stmEnergyCalib_h; - // MWD analysis variables + // PH analysis variables std::vector ADCs; // input waveform ADCs unsigned long int nADCs = 0; // size of input data unsigned long int i = 0; // iterator @@ -114,7 +115,7 @@ namespace mu2e { std::vector averaged_data; double sum = 0.0; // variable part of avarege() int count = 0; // counter variable - int16_t mwd_energy = 0; + int16_t ph_energy = 0; // Peak finding variables double baseline_mean = 0.0; @@ -145,7 +146,7 @@ namespace mu2e { STMMovingWindowDeconvolution::STMMovingWindowDeconvolution(const Parameters& conf) : art::EDProducer{conf}, _stmWaveformDigisToken(consumes(conf().stmWaveformDigisTag())), - channel(STMUtils::getChannel(conf().stmWaveformDigisTag())), + channel(STMChannel::findByName("HPGe")), // FIXME: don't hardcode this probably don't want to do what we had before and try to infer it from the art::InputTag like this "STMUtils::getChannel(config().stmWaveformDigisTag()))" tau(conf().tau()), M(conf().M()), L(conf().L()), @@ -153,17 +154,18 @@ namespace mu2e { thresholdgrad(conf().thresholdgrad()), defaultBaselineMean(conf().defaultBaselineMean()), defaultBaselineSD(conf().defaultBaselineSD()) { - produces(); + produces(); if (M < L) throw cet::exception("Configuration", "L (" + std::to_string(L) + ") is greater than M (" + std::to_string(M) + "), reconfigure\n"); verbosityLevel = conf().verbosityLevel() ? *(conf().verbosityLevel()) : 0; - if (verbosityLevel > 10) - verbosityLevel = 10; - _xAxis = conf().xAxis() ? *(conf().xAxis()) : ""; - makeTTreeMWD = conf().makeTTreeMWD() ? *(conf().makeTTreeMWD()) : false; + if (verbosityLevel >= 7){ + _xAxis = conf().xAxis() ? *(conf().xAxis()) : ""; + std::cout<<"BG"<<_xAxis< tfs; ttree = tfs->make("ttree", "STMMovingWindowDeconvolution pulse finding ttree"); ttree->Branch("ADC", &ADC, "ADC/S"); @@ -182,7 +184,7 @@ namespace mu2e { ttree->Branch("E", &E, "E/D"); ttree->Branch("waveformID", &waveformID, "waveformID/i"); }; - if (_xAxis != "") { + if (_xAxis == "") { if (verbosityLevel >= 5) { throw cet::exception("STMMovingWindowDecomposition") << "No xAxis scale defined despite requesting verbosity level >= 5" << std::endl; }; @@ -207,8 +209,9 @@ namespace mu2e { void STMMovingWindowDeconvolution::produce(art::Event& event) { // create output - std::unique_ptr outputMWDDigis(new STMMWDDigiCollection); + std::unique_ptr outputPHDigis(new STMPHDigiCollection); auto waveformDigisHandle = event.getValidHandle(_stmWaveformDigisToken); + // get prodition STMEnergyCalib const& stmEnergyCalib = _stmEnergyCalib_h.get(event.id()); pedestal = stmEnergyCalib.pedestal(channel); @@ -219,6 +222,13 @@ namespace mu2e { eventId = event.id().event(); waveformID = 0; for (STMWaveformDigi waveform : *waveformDigisHandle) { + + //Check prints -------------------------------------- + if (verbosityLevel>=7){ + std::cout << "event id = "<< event.id().event()<(peak_heights[i]); // When saturating the int16_t limit, deconvolution goes below the int16_t limit so the energy turns negative. This clips the energy and the limit - STMMWDDigi mwd_digi(peak_times[i], -1 * mwd_energy); // peak_heights are negative, make them positive here - if (mwd_digi.energy() < -100) + ph_energy = (peak_heights[i] < ADCMax) ? ADCMax : static_cast(peak_heights[i]); // When saturating the int16_t limit, deconvolution goes below the int16_t limit so the energy turns negative. This clips the energy and the limit + STMPHDigi ph_digi(peak_times[i], -1 * ph_energy); // peak_heights are negative, make them positive here + if (ph_digi.energy() < -100) throw cet::exception("logicError", "The peak height must be positive!"); - outputMWDDigis->push_back(mwd_digi); + outputPHDigis->push_back(ph_digi); if (makeTTreeEnergies) { - time = mwd_digi.time(); - E = mwd_digi.energy() * TTreeEnergyCalib; + time = ph_digi.time(); + E = ph_digi.energy() * TTreeEnergyCalib; ttree->Fill(); }; if (verbosityLevel > 3) - std::cout << "energy: " << mwd_digi.energy() << std::endl; + std::cout << "energy: " << ph_digi.energy() << std::endl; }; // Save data to TTree - if (makeTTreeMWD) { + if (makeTTreePH) { time = waveform.trigTimeOffset(); for (i = 0; i < nADCs; i++) { ADC = ADCs[i]; @@ -281,8 +291,8 @@ namespace mu2e { }; if (verbosityLevel) - std::cout << "MWD: " << channel.name() << ": " << outputMWDDigis->size() << " MWD digis found" << std::endl; - event.put(std::move(outputMWDDigis)); + std::cout << "MWD: " << channel.name() << ": " << outputPHDigis->size() << " PH digis found" << std::endl; + event.put(std::move(outputPHDigis)); }; @@ -292,6 +302,7 @@ namespace mu2e { for (int16_t data : ADCs) std::cout << data << ", "; std::cout << "\n" << std::endl; + // std::cout << "waveformhandle = " << waveformsDigisHandle<make(("h_baseline_mean_minus_stddev"+histsuffix.str()).c_str(), "Baseline Mean - StdDev", binning.nbins(),binning.low(),binning.high()); TH1D* h_peak_threshold = tfs->make(("h_peak_threshold"+histsuffix.str()).c_str(), "Threshold", binning.nbins(),binning.low(),binning.high()); + // now label them + h_waveform->GetXaxis()->SetTitle("Sample Number"); + h_waveform->GetYaxis()->SetTitle("ADCs"); + + h_deconvolved->GetXaxis()->SetTitle("Sample Number"); + h_deconvolved->GetYaxis()->SetTitle("ADCs"); + + h_differentiated->GetXaxis()->SetTitle("Sample Number"); + h_differentiated->GetYaxis()->SetTitle("ADCs"); + + h_averaged->GetXaxis()->SetTitle("Sample Number"); + h_averaged->GetYaxis()->SetTitle("ADCs"); + + h_baseline_mean->GetXaxis()->SetTitle("Sample Number"); + h_baseline_mean->GetYaxis()->SetTitle("ADCs"); + + h_baseline_mean_plus_stddev->GetXaxis()->SetTitle("Sample Number"); + h_baseline_mean_plus_stddev->GetYaxis()->SetTitle("ADCs"); + + h_baseline_mean_minus_stddev->GetXaxis()->SetTitle("Sample Number"); + h_baseline_mean_minus_stddev->GetYaxis()->SetTitle("ADCs"); + + h_peak_threshold->GetXaxis()->SetTitle("Sample Number"); + h_peak_threshold->GetYaxis()->SetTitle("ADCs"); + for (size_t i = 0; i < deconvolved_data.size(); ++i) { h_waveform->SetBinContent(i+1, waveform.adcs()[i] - pedestal); // remove the pedestal h_deconvolved->SetBinContent(i+1, deconvolved_data[i]); @@ -423,7 +459,12 @@ namespace mu2e { h_baseline_mean_minus_stddev->SetBinContent(i+1, baseline_mean - baseline_stddev); h_peak_threshold->SetBinContent(i+1, baseline_mean - nsigma_cut * baseline_stddev); } + TH1D* h_peaks = tfs->make(("h_peaks"+histsuffix.str()).c_str(), "Peaks", binning.nbins(),binning.low(),binning.high()); + + h_peaks->GetXaxis()->SetTitle("Sample Number"); + h_peaks->GetYaxis()->SetTitle("ADC peak"); + for (size_t i_peak = 0; i_peak < peak_heights.size(); ++i_peak) { // std::cout << "t = " << peak_times[i_peak] << ", E = " << peak_heights[i_peak] << std::endl; h_peaks->SetBinContent(peak_times[i_peak]+1, peak_heights[i_peak]); diff --git a/STMReco/test/plotDummyDigis.C b/STMReco/test/plotDummyDigis.C new file mode 100644 index 0000000000..4e2accf708 --- /dev/null +++ b/STMReco/test/plotDummyDigis.C @@ -0,0 +1,77 @@ +// +// +void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { + + TFile* file = new TFile(filename.c_str(), "READ"); + + TCanvas* c1 = new TCanvas(); + c1->Divide(2, 2); + + c1->cd(1); + TH1F* hRawWaveform = (TH1F*) file->Get("plotRawWaveformDigis/evt0_waveform0"); + TString title(hRawWaveform->GetTitle()); + title += " RAW"; + hRawWaveform->SetTitle(title); + hRawWaveform->SetLineWidth(3); + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->Draw("HIST"); + + c1->cd(2); + TH1F* hZSWaveform = (TH1F*) file->Get("plotZSWaveformDigis/evt0_waveform0"); + title = hZSWaveform->GetTitle(); + title += " ZS"; + hZSWaveform->SetTitle(title); + hZSWaveform->SetLineWidth(3); + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->Draw("HIST"); + + c1->cd(3); + TH1F* hPHWaveform = (TH1F*) file->Get("plotPHWaveformDigis/evt0_waveform0"); + title = hPHWaveform->GetTitle(); + title += " PH "; + hPHWaveform->SetTitle(title); + hPHWaveform->SetLineWidth(3); + hPHWaveform->SetLineColor(kRed); + hPHWaveform->Draw("HIST"); + +//New Canvas with all the histograms in the same axis +//No need to redefine the histograms from before + TCanvas* c2 = new TCanvas(); + c2->SetTitle("STMWaveformDigis Super Imposed"); + + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->SetLineWidth(3); + + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->SetLineWidth(3); + + hPHWaveform->SetLineColor(kRed); + hPHWaveform->SetLineWidth(3); + + hRawWaveform->GetYaxis()->SetRangeUser(-1800,1800); + hRawWaveform->Draw(); + hZSWaveform->Draw("SAME"); + hPHWaveform->Draw("SAME"); + + auto legend = new TLegend(); + legend->AddEntry(hRawWaveform, "Raw","1"); + legend->AddEntry(hZSWaveform,"ZS","1"); + legend->AddEntry(hPHWaveform,"PH","1"); + //legend->Draw(); + +//New Canvas where the histogram plots the number of bins from the other two + TCanvas* c3 = new TCanvas(); + c3->SetTitle("Number of bins from hPHWaveform and hZSWaveform"); + + TH1F* hDigiSize = new TH1F("hDigiSize","Number of bins in histograms",1000,0,1000); + hDigiSize->Fill(hPHWaveform->GetNbinsX()); + hDigiSize->Fill(hZSWaveform->GetNbinsX()); + hDigiSize->Draw(); + + + //Below is used to create a copy of the previous canvas + // TCanvas* c4 = new TCanvas(); + // c4->DrawClonePad(); + +} + diff --git a/STMReco/test/plotMWD.C b/STMReco/test/plotPH.c similarity index 78% rename from STMReco/test/plotMWD.C rename to STMReco/test/plotPH.c index 78c87b4473..849d44c477 100644 --- a/STMReco/test/plotMWD.C +++ b/STMReco/test/plotPH.c @@ -3,7 +3,7 @@ // This ROOT macro runs on the output of the STMMovingWindowDeconvolution module // with verbosity level set to 5 // -void plotMWD(std::string filename = "mwd.root") { +void plotPH(std::string filename = "ph.root") { TFile* file = new TFile(filename.c_str(), "READ"); @@ -13,42 +13,42 @@ void plotMWD(std::string filename = "mwd.root") { std::string evt_wvf = "evt1_wvf161"; - TH1F* hRawWaveform = (TH1F*) file->Get(("mwdLaBr/h_waveform_"+evt_wvf).c_str()); + TH1F* hRawWaveform = (TH1F*) file->Get(("phLaBr/h_waveform_"+evt_wvf).c_str()); hRawWaveform->SetStats(false); hRawWaveform->SetLineColor(kBlack); hRawWaveform->Draw("HIST"); leg->AddEntry(hRawWaveform, "raw waveform", "l"); - TH1F* hDeconvolved = (TH1F*) file->Get(("mwdLaBr/h_deconvolved_"+evt_wvf).c_str()); + TH1F* hDeconvolved = (TH1F*) file->Get(("phLaBr/h_deconvolved_"+evt_wvf).c_str()); hDeconvolved->SetLineColor(kOrange); hDeconvolved->Draw("HIST SAME"); leg->AddEntry(hDeconvolved, "deconvolution (tau parameter)", "l"); - TH1F* hDifferentiated = (TH1F*) file->Get(("mwdLaBr/h_differentiated_"+evt_wvf).c_str()); + TH1F* hDifferentiated = (TH1F*) file->Get(("phLaBr/h_differentiated_"+evt_wvf).c_str()); hDifferentiated->SetLineColor(kSpring); hDifferentiated->Draw("HIST SAME"); leg->AddEntry(hDifferentiated, "differentiation (M parameter)", "l"); - TH1F* hAveraged = (TH1F*) file->Get(("mwdLaBr/h_averaged_"+evt_wvf).c_str()); + TH1F* hAveraged = (TH1F*) file->Get(("phLaBr/h_averaged_"+evt_wvf).c_str()); hAveraged->SetLineColor(kRed); hAveraged->SetLineWidth(2); hAveraged->Draw("HIST SAME"); leg->AddEntry(hAveraged, "average (L parameter)", "l"); - TH1F* hBaselineMean = (TH1F*) file->Get(("mwdLaBr/h_baseline_mean_"+evt_wvf).c_str()); + TH1F* hBaselineMean = (TH1F*) file->Get(("phLaBr/h_baseline_mean_"+evt_wvf).c_str()); hBaselineMean->SetLineColor(kBlue); hBaselineMean->SetLineWidth(2); hBaselineMean->Draw("HIST SAME"); leg->AddEntry(hBaselineMean, "baseline mean", "l"); - TH1F* hPeaks = (TH1F*) file->Get(("mwdLaBr/h_peaks_"+evt_wvf).c_str()); + TH1F* hPeaks = (TH1F*) file->Get(("phLaBr/h_peaks_"+evt_wvf).c_str()); hPeaks->SetLineColor(kBlue); hPeaks->SetLineWidth(2); hPeaks->Draw("HIST SAME"); leg->AddEntry(hPeaks, "found peaks", "l"); - TH1F* hPeakThresh = (TH1F*) file->Get(("mwdLaBr/h_peak_threshold_"+evt_wvf).c_str()); + TH1F* hPeakThresh = (TH1F*) file->Get(("phLaBr/h_peak_threshold_"+evt_wvf).c_str()); hPeakThresh->SetLineColor(kRed); hPeakThresh->SetLineWidth(2); hPeakThresh->SetLineStyle(kDashed);