diff --git a/Analyses/src/CountVirtualDetectorHits_module.cc b/Analyses/src/CountVirtualDetectorHits_module.cc index 656f5db9a3..2f80b510c6 100644 --- a/Analyses/src/CountVirtualDetectorHits_module.cc +++ b/Analyses/src/CountVirtualDetectorHits_module.cc @@ -51,10 +51,17 @@ namespace mu2e { : art::EDAnalyzer(conf), StepPointMCsToken(consumes(conf().stepPointMCsTag())), enabledVDs(conf().enabledVDs()) { - // Create list of unique enabled virtual detectors - std::set enabledVDsSet(enabledVDs.begin(), enabledVDs.end()); + // Create list of unique enabled virtual detectors, preserving the order + std::vector unqiueEnabledVDsVec;; + std::unordered_set uniqueEnabledVDsSet; + for (const int & enabledVD : enabledVDs) { + if (uniqueEnabledVDsSet.find(enabledVD) == uniqueEnabledVDsSet.end()) { + uniqueEnabledVDsSet.insert(enabledVD); + unqiueEnabledVDsVec.push_back(enabledVD); + }; + }; enabledVDs.clear(); - enabledVDs.insert(enabledVDs.end(), enabledVDsSet.begin(), enabledVDsSet.end()); + enabledVDs = unqiueEnabledVDsVec; // Insert _enabledVDs.size() zeros into the counter vector nVDs = enabledVDs.size(); diff --git a/BTrkLegacy/inc/TrkErrCode.hh b/BTrkLegacy/inc/TrkErrCode.hh index 48ed822520..80062b2868 100644 --- a/BTrkLegacy/inc/TrkErrCode.hh +++ b/BTrkLegacy/inc/TrkErrCode.hh @@ -98,5 +98,3 @@ private: std::ostream& operator<<(std::ostream& os, const TrkErrCode& trkerr); #endif - - diff --git a/BTrkLegacy/inc/TrkT0.hh b/BTrkLegacy/inc/TrkT0.hh index 621f3d1bc1..9caeb2ed2f 100644 --- a/BTrkLegacy/inc/TrkT0.hh +++ b/BTrkLegacy/inc/TrkT0.hh @@ -13,4 +13,4 @@ struct TrkT0 { double _t0err; // estimate of t0 error }; -#endif \ No newline at end of file +#endif diff --git a/DAQ/src/MSDHitsFromDTCEvents_module.cc b/DAQ/src/MSDHitsFromDTCEvents_module.cc index 1e887a0895..3a53e38246 100644 --- a/DAQ/src/MSDHitsFromDTCEvents_module.cc +++ b/DAQ/src/MSDHitsFromDTCEvents_module.cc @@ -193,4 +193,3 @@ artdaq::Fragments art::MSDHitsFromDTCEvents::getFragments(art::Event& event) { DEFINE_ART_MODULE(art::MSDHitsFromDTCEvents) // ====================================================================== - diff --git a/DAQ/src/STMBinaryDigisFromFragments_module.cc b/DAQ/src/STMBinaryDigisFromFragments_module.cc index d3c0bf4304..d3c0f0c35b 100644 --- a/DAQ/src/STMBinaryDigisFromFragments_module.cc +++ b/DAQ/src/STMBinaryDigisFromFragments_module.cc @@ -41,11 +41,11 @@ class art::STMBinaryDigisFromFragments : public EDProducer 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 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 @@ -101,7 +101,7 @@ STMBinaryDigisFromFragments::STMBinaryDigisFromFragments(const art::EDProducer:: if(!_rawOut || !_zsOut || !_phOut || !_rawHeaderOut || !_eventOut){ throw cet::exception("FILEOPEN")<< "Failed to open one or more output files\n"; } - + } STMBinaryDigisFromFragments::~STMBinaryDigisFromFragments(){ @@ -114,16 +114,16 @@ STMBinaryDigisFromFragments::~STMBinaryDigisFromFragments(){ if( _rawHeaderOut.is_open() ) _rawHeaderOut.close(); if( _eventOut.is_open() ) - _eventOut.close(); + _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); @@ -145,25 +145,25 @@ void STMBinaryDigisFromFragments::produce(Event& event) size_t localRaw_frags{0}; size_t localZS_frags{0}; size_t localPH_frags{0}; - size_t zeroRaw_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){ @@ -173,7 +173,7 @@ void STMBinaryDigisFromFragments::produce(Event& event) _totalInner += blocks; //loop over container - // i index corresponds to inner frag + // i index corresponds to inner frag for (size_t i = 0; i < cont_frag.block_count(); ++i){ auto inner_frag = cont_frag.at(i); @@ -185,19 +185,19 @@ void STMBinaryDigisFromFragments::produce(Event& event) if (stm_frag.isRaw()) { ++_totalRaw; //Increment job counter - ++localRaw_frags;//Increment event 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 + 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){ @@ -209,172 +209,172 @@ void STMBinaryDigisFromFragments::produce(Event& event) //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 + 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"; - } - } + 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); - + //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; + 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"; - // } - + ++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(); + ++_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; + writePayload(_phOut,digiPtr, digiWords); + extractedPHDigis += digiWords; - }//End of isPH and is checks + }//End of isPH and is checks //---Combined stream write w. order preserved ---- { @@ -412,7 +412,7 @@ void STMBinaryDigisFromFragments::produce(Event& event) } } } //End of frags loop---George suggestions - if (_verbosityLevel >= 2 ){ + 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"; @@ -423,7 +423,7 @@ void STMBinaryDigisFromFragments::produce(Event& event) 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"; @@ -432,10 +432,10 @@ void STMBinaryDigisFromFragments::produce(Event& event) 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() // ====================================================================== diff --git a/DAQ/src/STMPrintFragments_module.cc b/DAQ/src/STMPrintFragments_module.cc index adf5f8d0b1..9c909ed54f 100644 --- a/DAQ/src/STMPrintFragments_module.cc +++ b/DAQ/src/STMPrintFragments_module.cc @@ -69,11 +69,11 @@ void STMPrintFragments::analyze(const Event& event) 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(); @@ -83,9 +83,9 @@ void STMPrintFragments::analyze(const Event& event) const auto stmDataEnd = reinterpret_cast(dataEnd); std::cout << "Frag_ID = " << frag_id << std::endl; std::cout << "Container block_count = "< +#include // art includes #include "art/Framework/Core/EDProducer.h" @@ -40,7 +41,11 @@ namespace mu2e { fhicl::Atom z{ Name("z"), Comment("z position of generated photon")}; fhicl::OptionalAtom px{ Name("px"), Comment("x momentum of generated photon")}; fhicl::OptionalAtom py{ Name("py"), Comment("y momentum of generated photon")}; + fhicl::OptionalAtom deltax{ Name("deltax"), Comment("Difference in x position of generated photon, use to calculate targeted position")}; + fhicl::OptionalAtom deltay{ Name("deltay"), Comment("Difference in y position of generated photon, use to calculate targeted position")}; + fhicl::OptionalAtom deltaz{ Name("deltaz"), Comment("Difference in z position of generated photon, use to calculate targeted position")}; fhicl::Atom E{ Name("E"), Comment("Energy of generated photon")}; + fhicl::OptionalAtom verbose{ Name("verbose"), Comment("Print verbose messages")}; }; using Parameters = art::EDProducer::Table; explicit PhotonGun(const Parameters& conf); @@ -49,7 +54,9 @@ namespace mu2e { private: double x = 0.0, y = 0.0, z = 0.0; double px = 0.0, py = 0.0, pz = 0.0; + double deltax = 0.0, deltay = 0.0, deltaz = 0.0; double E = 0.0; + bool verbose = false; }; PhotonGun::PhotonGun(const Parameters& conf): @@ -59,18 +66,33 @@ namespace mu2e { z(conf().z()), E(conf().E()) { produces(); - produces(); + if (E < std::numeric_limits::epsilon()) + throw cet::exception("RANGE") << "Energy must be greater than zero, exiting."; px = conf().px() ? *conf().px() : 0; - px = conf().py() ? *conf().py() : 0; + py = conf().py() ? *conf().py() : 0; if ((px*px + py*py) > (E*E)) throw cet::exception("RANGE") << "magnitude of px and py is greater than E, exiting."; pz = std::sqrt(E*E - px*px - py*py); + deltax = conf().deltax() ? *conf().deltax() : 0; + deltay = conf().deltay() ? *conf().deltay() : 0; + deltaz = conf().deltaz() ? *conf().deltaz() : 0; + if ((std::abs(px) > std::numeric_limits::epsilon() || std::abs(py) > std::numeric_limits::epsilon()) && (std::abs(deltax) > std::numeric_limits::epsilon() || std::abs(deltay) > std::numeric_limits::epsilon() || std::abs(deltaz) > std::numeric_limits::epsilon())) + throw cet::exception("RANGE") << "Cannot specify both momentum and delta position, exiting."; + verbose = conf().verbose() ? *conf().verbose() : false; }; void PhotonGun::produce(art::Event& event) { std::unique_ptr output(new GenParticleCollection); const CLHEP::Hep3Vector pos(x, y, z); - const CLHEP::Hep3Vector p(px, py, pz); + CLHEP::Hep3Vector p(px, py, pz); + if (std::abs(deltax) > std::numeric_limits::epsilon() || + std::abs(deltay) > std::numeric_limits::epsilon() || + std::abs(deltaz) > std::numeric_limits::epsilon()) { + if (verbose) + std::cout << "PhotonGun: Using delta position to calculate momentum." << std::endl; + const CLHEP::Hep3Vector delta(deltax, deltay, deltaz); + p = delta.unit() * E; + } CLHEP::HepLorentzVector mom(p, E); output->push_back(GenParticle(PDGCode::gamma, GenId::particleGun, pos, mom, 0.)); event.put(std::move(output)); diff --git a/RecoDataProducts/inc/MSDHit.hh b/RecoDataProducts/inc/MSDHit.hh index ba4d8707fd..fdab4b1f82 100644 --- a/RecoDataProducts/inc/MSDHit.hh +++ b/RecoDataProducts/inc/MSDHit.hh @@ -33,4 +33,3 @@ typedef std::vector MSDHitCollection; } // namespace mu2e #endif /* RecoDataProducts_MSDHit_hh */ - diff --git a/RecoDataProducts/inc/STMWaveformDigi.hh b/RecoDataProducts/inc/STMWaveformDigi.hh index 6d0dc8551c..97a16a05aa 100644 --- a/RecoDataProducts/inc/STMWaveformDigi.hh +++ b/RecoDataProducts/inc/STMWaveformDigi.hh @@ -18,12 +18,12 @@ namespace mu2e { class STMWaveformDigi { public: - // Initialise all variables + // Initialise all variables STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _adcs(std::vector()) {}; - // Constructor for timing plus trig offset + // 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 + // 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; } diff --git a/STMMC/CMakeLists.txt b/STMMC/CMakeLists.txt index f656f6d24a..a3bb1bdb02 100644 --- a/STMMC/CMakeLists.txt +++ b/STMMC/CMakeLists.txt @@ -57,6 +57,22 @@ cet_build_plugin(ShiftVirtualDetectorStepPointMCs art::module Offline::MCDataProducts ) +cet_build_plugin(SimParticleDump art::module + REG_SOURCE src/SimParticleDump_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::GlobalConstantsService + Offline::MCDataProducts +) + +cet_build_plugin(SimParticleAndVDBacktrace art::module + REG_SOURCE src/SimParticleAndVDBacktrace_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::GlobalConstantsService + Offline::MCDataProducts +) + cet_build_plugin(STMResamplingProducer art::module REG_SOURCE src/STMResamplingProducer_module.cc LIBRARIES REG diff --git a/STMMC/fcl/Absorber.fcl b/STMMC/fcl/Absorber.fcl index 14c1785f8d..bf07469418 100644 --- a/STMMC/fcl/Absorber.fcl +++ b/STMMC/fcl/Absorber.fcl @@ -25,18 +25,18 @@ physics: { producers : { generate : { module_type : PhotonGun - x : @local::Efficiency.Absorber.x - y : @local::Efficiency.Absorber.y - z : @local::Efficiency.Absorber.z + x : @local::ComponentPositions.Absorber.x + y : @local::ComponentPositions.Absorber.y + z : @local::ComponentPositions.Absorber.z E : @local::Efficiency.PhotonEnergy + verbose : false } g4run : @local::g4run } analyzers : { EDep : { - module_type : MakeVirtualDetectorTree - VirtualDetectorId : 101 + module_type : VirtualDetectorTree StepPointMCsTag: "g4run:virtualdetector" SimParticlemvTag: "g4run:" } @@ -50,7 +50,7 @@ physics: { physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] -services.SeedService.baseSeed : 8 +services.SeedService.baseSeed : 18 services.SeedService.maxUniqueEngines : 20 services.TFileService.fileName : @local::Efficiency.OutputFilename # physics.producers.g4run.debug.trackingVerbosityLevel : 1 diff --git a/STMMC/fcl/AbsorberFromSTHPGe.fcl b/STMMC/fcl/AbsorberFromSTHPGe.fcl new file mode 100644 index 0000000000..db1cfd9f34 --- /dev/null +++ b/STMMC/fcl/AbsorberFromSTHPGe.fcl @@ -0,0 +1,60 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates photons being fired from the ST incident on the STM detectors to investigate the relative shift in spectrum +# Original author : Pawel Plesniak +# Note - edit the parameters in this prolog.fcl, aside from which detector is being used which is stored in the generate module + +process_name: HPGeAbsorberSpectrumShift +source : { + module_type : EmptyEvent + maxEvents : @local::Efficiency.NPhotons +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::Efficiency.ST.x + y : @local::Efficiency.ST.y + z : @local::Efficiency.ST.z + deltax : @local::Efficiency.FromSTToDet.HPGe.x + deltay : @local::Efficiency.FromSTToDet.HPGe.y + deltaz : @local::Efficiency.FromSTToDet.HPGe.z + E : @local::Efficiency.PhotonEnergy + verbose: false + } + g4run : @local::g4run + } + + analyzers : { + EDep : { + module_type : VirtualDetectorTree + StepPointMCsTag: "g4run:virtualdetector" + SimParticlemvTag: "g4run:" + } + } + + p1 : [ generate, g4run ] + trigger_paths : [p1] + o1 : [ EDep ] + end_paths: [ o1 ] +} + +physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] +services.SeedService.baseSeed : 8 +services.SeedService.maxUniqueEngines : 20 +services.TFileService.fileName : @local::Efficiency.OutputFilename +# physics.producers.g4run.debug.trackingVerbosityLevel : 1 +# physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/AbsorberFromSTLaBr.fcl b/STMMC/fcl/AbsorberFromSTLaBr.fcl new file mode 100644 index 0000000000..a128c021ce --- /dev/null +++ b/STMMC/fcl/AbsorberFromSTLaBr.fcl @@ -0,0 +1,60 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates photons being fired from the ST incident on the STM detectors to investigate the relative shift in spectrum +# Original author : Pawel Plesniak +# Note - edit the parameters in this prolog.fcl, aside from which detector is being used which is stored in the generate module + +process_name: HPGeAbsorberSpectrumShift +source : { + module_type : EmptyEvent + maxEvents : @local::Efficiency.NPhotons +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::Efficiency.ST.x + y : @local::Efficiency.ST.y + z : @local::Efficiency.ST.z + deltax : @local::Efficiency.FromSTToDet.LaBr.x + deltay : @local::Efficiency.FromSTToDet.LaBr.y + deltaz : @local::Efficiency.FromSTToDet.LaBr.z + E : @local::Efficiency.PhotonEnergy + verbose: false + } + g4run : @local::g4run + } + + analyzers : { + EDep : { + module_type : VirtualDetectorTree + StepPointMCsTag: "g4run:virtualdetector" + SimParticlemvTag: "g4run:" + } + } + + p1 : [ generate, g4run ] + trigger_paths : [p1] + o1 : [ EDep ] + end_paths: [ o1 ] +} + +physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] +services.SeedService.baseSeed : 8 +services.SeedService.maxUniqueEngines : 20 +services.TFileService.fileName : @local::Efficiency.OutputFilename +# physics.producers.g4run.debug.trackingVerbosityLevel : 1 +# physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/EventFilter.fcl b/STMMC/fcl/EventFilter.fcl new file mode 100644 index 0000000000..e17ab65e2b --- /dev/null +++ b/STMMC/fcl/EventFilter.fcl @@ -0,0 +1,42 @@ +# Script to generate a separate file associated with a single event, for testing purposes. This is not intended to be used in production, but it can be useful for debugging and development. +# The event ID is specified in the 'idsToMatch' parameter of the EventIDFilter module. In this example, it is set to "1:0:5", which corresponds to run 1, subrun 0, event 5. You can change this value to select a different event. +# Mainly for use with the output of FirePhotonToDetector.fcl, as input to HPGeWaveformGenerationAndAnalysis.fcl, but it can be used with any input file that contains the specified event. +# Original author: Pawel Plesniak + +#include "Offline/fcl/minimalMessageService.fcl" + +process_name: STMCat + +source: { + module_type: RootInput +} + +physics: { + filters: { + isolateEvent: { + module_type: EventIDFilter + # The framework specifically asked for 'idsToMatch' + idsToMatch: [ "1:0:5" ] + } + } + + # This is the path name that SelectEvents must reference + STMCompressedPath: [ isolateEvent ] + trigger_paths: [ STMCompressedPath ] + o1 : [ out ] + end_paths: [ o1 ] +} + +services: { + message: @local::default_message +} + +outputs: { + out: { + module_type: RootOutput + fileName: "Filtered.art" + SelectEvents: ["STMCompressedPath"] + } +} + +services.scheduler.wantSummary: true diff --git a/STMMC/fcl/FirePhotonToDetector.fcl b/STMMC/fcl/FirePhotonToDetector.fcl new file mode 100644 index 0000000000..19e968ca10 --- /dev/null +++ b/STMMC/fcl/FirePhotonToDetector.fcl @@ -0,0 +1,60 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates N photon(s) being fired upstream of the HPGe detector to validate the behaviour of the digitizer and analysis tools +# Original author : Pawel Plesniak + +process_name: FirePhotonToDetector +source : { + module_type : EmptyEvent + maxEvents : 1 +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::DigitizationTester.HPGe.x + y : @local::DigitizationTester.HPGe.y + z : @local::DigitizationTester.HPGe.z + E : 1.809 #@nil + verbose : true + } + g4run : @local::g4run + } + + p1 : [ generate, g4run ] + trigger_paths : [ p1 ] + o1 : [ out ] + end_paths: [ o1 ] +} + +outputs : { + out : { + module_type : RootOutput + fileName : "FirePhotonToDetector.art" + outputCommands: [ + "drop *_*_*_*", + "keep mu2e::StepPointMCs_g4run_STMDet_FirePhotonToDetector", + "keep mu2e::SimParticlemv_g4run__FirePhotonToDetector" + ] + } +} + +# physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet] +services.SeedService.baseSeed : 18 +services.SeedService.maxUniqueEngines : 20 +# services.TFileService.fileName : @local::Efficiency.OutputFilename +physics.producers.g4run.debug.trackingVerbosityLevel : 1 +physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/HPGeReco.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl similarity index 67% rename from STMMC/fcl/HPGeReco.fcl rename to STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index 4627443482..d3034112a2 100644 --- a/STMMC/fcl/HPGeReco.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -1,4 +1,13 @@ - +# Executes the full HPGe waveform generation and analysis chain on STMMC data +# Outline: +# 1. DigiHPGe - generates HPGe waveforms from StepPointMCs +# 2. concatenateWaveformsHPGe - concatenates microspill waveforms into spill waveforms +# 3. (optional) ZSHPGe - zero-suppresses the concatenated waveforms +# 4. MWDHPGe - performs moving window deconvolution on the concatenated (or zero-suppressed) waveforms +# 5. Ou1tputs the MWD digis to a ROOT file +# Note - each module is documented in the _module.cc file. This is the clearest way to understand what each module does and what parameters it takes. +# Usage: +# 1. Select which data you want to run with or without zero suppression by replacing "@nil" in trigger_paths (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. trigger_paths: [digitization_path_ZS] # Original author: Pawel Plesniak #include "Offline/fcl/standardServices.fcl" @@ -17,6 +26,7 @@ services : @local::Services.Sim physics : { producers : { DigiHPGe : { + # Digitize the StepPointMCs from the HPGe crystals into waveforms module_type : HPGeWaveformsFromStepPointMCs StepPointMCsTagEle : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagEle StepPointMCsTagMu : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagMu @@ -27,17 +37,20 @@ physics : { risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount makeTTree : false + makeADCPlot : false # Note - this does not work, I don't have the time to fix this. timeOffset : 0 resetEventNumber : @local::HPGeDigitization.resetEventNumber - verbosityLevel : 0 + verbosityLevel : 1 } concatenateWaveformsHPGe : { + # Concatenate the microspill waveforms into macrospill waveforms module_type : ConcatenateDigitizedWaveforms STMWaveformDigisTag : @local::HPGeDigitization.concatenation.STMWaveformDigisTag nMerge : @local::HPGeDigitization.concatenation.nMerge makeTTree : false } ZSHPGe : { + # Apply zero-suppression to the concatenated waveforms module_type : STMZeroSuppression stmWaveformDigisTag : @local::STMMCAnalysis.ZS.HPGe.stmWaveformDigisTag.concatenated tbefore : @local::STMMCAnalysis.ZS.HPGe.tbefore @@ -45,11 +58,12 @@ physics : { threshold : @local::STMMCAnalysis.ZS.HPGe.threshold window : @local::STMMCAnalysis.ZS.HPGe.window naverage : @local::STMMCAnalysis.ZS.HPGe.naverage - verbosityLevel : 0 + verbosityLevel : 1 makeTTreeGradients: false makeTTreeWaveforms: false } - PHHPGe : { + MWDHPGe : { + # Perform moving window deconvolution on the concatenated (or zero-suppressed) waveforms to determine the pulse energies module_type : STMMovingWindowDeconvolution stmWaveformDigisTag : @local::STMMCAnalysis.PH.HPGe.stmWaveformDigisTag.concatenated tau : @local::STMMCAnalysis.MWD.HPGe.tau @@ -62,7 +76,7 @@ physics : { makeTTreePH: false makeTTreeEnergies: false TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin - verbosityLevel : 0 + verbosityLevel : 1 xAxis : "" } } @@ -72,8 +86,9 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } - digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, PHHPGe] - trigger_paths : [digitization_path] + digitization_path_no_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, MWDHPGe] + digitization_path_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] + trigger_paths : @nil # Populate me! output_path : [compressedOutput] end_paths : [output_path] } @@ -82,7 +97,7 @@ outputs : { compressedOutput : { module_type : RootOutput fileName : "dts.owner.HPGeReco.version.sequencer.art" - SelectEvents: ["digitization_path"] + SelectEvents: @nil # Populate me! outputCommands: [ "drop *_*_*_*", "keep mu2e::STMPHDigis_PHHPGe_*_HPGeReco" diff --git a/STMMC/fcl/ROOTAnalysisDump.fcl b/STMMC/fcl/ROOTAnalysisDump.fcl new file mode 100644 index 0000000000..aa899682ef --- /dev/null +++ b/STMMC/fcl/ROOTAnalysisDump.fcl @@ -0,0 +1,107 @@ +# Extracts parameters from STM simulation scripts into ROOT TTrees for use in the scripts defined it Mu2e/STMAnalysis repository +# This generates ROOT files containing different things depending on which analyzer is selected in the o1 path: +# - stage1virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 1 simulation +# - stage2virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 2 simulation +# - HPGeWaveforms : Dumps HPGe waveforms from Stage 2 simulation +# - LaBrWaveforms : Dumps LaBr3 waveforms from Stage 2 simulation (not implemented yet) +# - ZSwaveforms : Dumps zero-suppressed waveforms from HPGe digitisers in Stage 2 simulation (not implemented yet) +# - MWDspectra : Dumps MWD spectra from HPGe digitisers in Stage 2 simulation +# Usage: +# 1. Select which data you want to dump to a ROOT TTree by replacing "@nil" in the list o1 (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. o1: [stage1virtualdetectorStepPointMCs] +# 2. Set the output ROOT file name in services.TFileService.fileName +# 3. Prepare a text file containing a list of input STMMC data files (one file name per line) +# 3. Run as: `mu2e -c ROOTAnalysisDump.fcl -S FileWithListOfDataFiles.txt` +# +# Original author: Pawel Plesniak + +#include "Offline/STMMC/fcl/prolog.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name: ROOTAnalysisDump + +source : { + module_type : RootInput + fileNames: @nil +} +services : { + message : @local::default_message + GlobalConstantsService : { + inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" + } +} +physics: { + analyzers : { + # Simulation stage 1 ROOT files + Stage1virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage1.StepPointMCsTag # 1809 + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage1SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + } + + # Simulation stage 2 ROOT files + Stage2virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagVirtualdetector + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage2SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + } + # TODO - see if this can be removed. + # ROOTDump : { + # module_type : SignalParticleBacktraceVDs + # # StepPointMCsTag : "compressDetStepMCsSTM:virtualdetector" #Stage 1 + # # # StepPointMCsTag: "compressDetStepMCsSTM:" #Stage 1 1809 as Stage 1 Mu3 + # # StartVirtualDetectorId : 101 + # # TrackVirtualDetectorIDs : [101, 100, 10, 9, 8, 2, 1] + # StepPointMCsTag : "compressSTMDet:virtualdetector" #Stage 2 + # StartVirtualDetectorId : 90 + # TrackVirtualDetectorIDs : [90, 100, 10, 9, 8, 2, 1] + # EnergyWindowSize : 0.1 + # } + + # HPGe energy depositions in crystal - MC truth + HPGeEDeps : { + module_type : HPGeTree + StepPointMCsTag : "compressSTMDet:STMDet" + SimParticlemvTag : "compressSTMDet:" + Detector: @local::ROOTAnalysisDump.Stage2.DetectorName.HPGe + } +# HPGeZSWaveforms : { +# module_type : @nil # Doesn't exist yet! +# StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet +# SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag # TODO: Remove the SimParticle dependency +# Detector: @local::SimplifyStage2Data.DetectorName.HPGe +# } + HPGeMWDspectra : { + module_type : MWDTree + STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag + EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin # Once the proditions are implemneted, this should be removed. + } + + # LaBr3 waveform ROOT files - not implemented yet + LaBrWaveforms : { + module_type : @nil + } + LaBrZSWaveforms : { + module_type : @nil + } + LaBrMWDspectra : { + module_type : @nil + } + } + + o1 : ["HPGeEDeps"] # Populate me! + end_paths: [o1] +} + +services.TFileService.fileName : "S2E.HPGeTree.root" # Populate me! diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 2dca0648fa..275b17d2d8 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -1,10 +1,72 @@ -# -# prolog.fcl for fcl paramaters that will be used in STM simulation modules -# +# Defines standard paramaters used in STM simulation modules +# Adapted by: Pawel Plesniak + +# TODOs before Pawel uploads: +# - Review all parameters for correctness +# - Add comments where necessary +# - Verify that all necessary parameters are included +# - Simplify the ROOTAnalysisDump data product names by defining them separately here, once the simulation generats them correctly +# - SimParticlemv -> SimParticles +# - StepPointMCsTag -> StepPointMCs #include "Offline/STMReco/fcl/prolog.fcl" +#include "Production/JobConfig/pileup/STM/prolog.fcl" BEGIN_PROLOG +# Data product definitions +DataProducts : { + # Data products generated in Stage 1 of the simulation + Stage1 : { + StepPointMCs : @local::STMSimDataProducts.Stage1.CompressedOutput.StepPointMCs + SimParticles : @local::STMSimDataProducts.Stage1.CompressedOutput.SimParticles + EndVirtualdetectorID : @local::STMPileup.ResamplingProducer.VirtualDetectorID + } + + # Data products generated in Stage 2 of the simulation + Stage2 : { + StepPointMCs : { + Virtualdetector : @local::STMSimDataProducts.Stage2.CompressedOutput.StepPointMCs.Virtualdetector + Detector : @local::STMSimDataProducts.Stage2.CompressedOutput.StepPointMCs.Detector + } + SimParticles : @local::STMSimDataProducts.Stage2.CompressedOutput.SimParticles + EndVirtualdetectorID : { + HPGe : 90 + LaBr : 89 + } + } + + # Data products used in mixing from POT frame to microspill frame + MixingStepPointMCTags : { + StepPointMCsTagEle : "STMStepMixerEle:STMDet" + StepPointMCsTagMu : "STMStepMixerMu:STMDet" + StepPointMCsTag1809 : "STMStepMixer1809:STMDet" + } + + # Data products generated in STM waveform generation and analysis + Waveforms : { + HPGe : { + uSpill : "DigiHPGe:" + concatenated : "concatenateWaveformsHPGe:" + ZS : "ZSHPGe:" + } + LaBr : { + uSpill : "DigiLaBr:" + concatenated : "concatenateWaveformsLaBr:" + ZS : "ZSLaBr:" + } + } + + # Data products generated in STM MWD analysis + PeakEnergies : { + HPGe : "MWDHPGe:" + LaBr : "MWDLaBr:" + } + + # Miscellaneous parameters + nMicrospillsPerSpill : 31858 +} + +## Simulation parameters # Geometry, all in mm ComponentPositions : { # In beam order Beam : { @@ -13,7 +75,7 @@ ComponentPositions : { # In beam order } Absorber : { # 50mm x 50mm x 390mm x : -3944.6 - y : -24.9 # TODO - why is this not 0.0? + y : 0 z : 39900 } VD101 : { # Centre @@ -34,76 +96,39 @@ ComponentPositions : { # In beam order LaBr : { # aperture centre x : -3863.4 y : 0.0 + z : 40404 } -} - -# Stage 1 propagation parameters -ResamplingProducer : { - StepPointMCsTag : "g4run:virtualdetector" - VirtualDetectorID : 101 -} -ResamplingFilter : { - StepPointMCsTag : "compressDetStepMCsSTM" -} -VirtualDetectorCounter : { - StepPointMCsTag : "g4run:virtualdetector" - virtualDetectorIDs : [88, 89, 90, 100, 101] -} - -# ShiftVirtualDetectorStepPointMCs_module.cc parameters -ShiftVD101Steps : { - StepPointMCsTag : "compressDetStepMCsSTM:virtualdetector" - VirtualDetectorID : 101 - InputRadius : 200 # VD101 aperture radius - OutputRadius : 3.98942280401 # SSC aperture radius - HPGeUpStr : { # Position just upstream of absorber centered on the HPGe SSC aperture - x : -3944.6 - y : 0.0 - z : 39906.0 - } - LaBrUpStr : { # Position just upstream of SSC centered on the Labr SSC aperture - x : -3863.4 + ST : { # ST centre + x : -3904.0 y : 0.0 - z : 40300.9 + z : 627.0 } - pdgID : 0 } -# Mixing parameters -MixSTMEvents : { - extendedMean2BB : 3.93e7 # Copied from Production/JobConfig/mixing/TwoBB.fcl - cutMax2BB : 2.36e8 # Copied from Production/JobConfig/mixing/TwoBB.fcl - extendedMean1BB : 1.58e7 # Copied from Production/JobConfig/mixing/OneBB.fcl - cutMax1BB : 9.48e7 # Copied from Production/JobConfig/mixing/OneBB.fcl - SDF : 0.6 - nPOTs : 26919873604557116 - nMicroSpills : 690253169 - meanEventsPerPOTFactors : { - EleBeamCat : 5.465876332164181e-11 - MuBeamCat : 7.887481312841522e-13 - TargetStopsCat1809 : 3.114418308224206e-16 - } -} +# Waveform generation and digitization parameters STMMCAnalysis : { MixedEventsTags : { - StepPointMCsTagEle : "STMStepMixerEle:STMDet" - StepPointMCsTagMu : "STMStepMixerMu:STMDet" - StepPointMCsTag1809 : "STMStepMixer1809:STMDet" + StepPointMCsTagEle : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTagEle + StepPointMCsTagMu : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTagMu + StepPointMCsTag1809 : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTag1809 } ZS : { HPGe : { stmWaveformDigisTag : { - concatenated : "concatenateWaveformsHPGe:" - uSpill : "DigiHPGe:" + uSpill : @local::DataProducts.Waveforms.HPGe.uSpill + concatenated : @local::DataProducts.Waveforms.HPGe.concatenated } tbefore : @local::STM.HPGe.tbefore tafter : @local::STM.HPGe.tafter threshold : @local::STM.HPGe.threshold - window : @local::STM.HPGe.window # TODO - when using the optimized number from STMReco, the gradients were fluctuating between the minimum band of an int16_t and 0, but this disappeared when set to 50. Figure out why. + window : @local::STM.HPGe.window naverage : @local::STM.HPGe.naverage } LaBr : { - stmWaveformDigisTag : "DigiLaBr:" + stmWaveformDigisTag : { + uSpill : @local::DataProducts.Waveforms.LaBr.uSpill + concatenated : @local::DataProducts.Waveforms.LaBr.concatenated + } tbefore : @local::STM.LaBr.tbefore tafter : @local::STM.LaBr.tafter threshold : @local::STM.HPGe.threshold @@ -114,8 +139,8 @@ STMMCAnalysis : { PH : { HPGe : { stmWaveformDigisTag : { - ZS: "ZSHPGe:" - concatenated: "concatenateWaveformsHPGe:" + concatenated: @local::DataProducts.Waveforms.HPGe.concatenated + ZS: @local::DataProducts.Waveforms.HPGe.ZS } tau : @local::STM.HPGe.tau M : @local::STM.HPGe.M @@ -130,40 +155,96 @@ STMMCAnalysis : { full: @local::STM.HPGe.defaultBaselineSD suppressed: 0 } - STMPHDigiTag : "PHHPGe:" + STMMWDDigiTag : @local::DataProducts.PeakEnergies.HPGe } LaBr : { - stmWaveformDigisTag : "ZSLaBr:" + stmWaveformDigisTag : { + concatenated: @local::DataProducts.Waveforms.LaBr.concatenated + ZS: @local::DataProducts.Waveforms.LaBr.ZS + } 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 - STMPHDigiTag : "PHLaBr:" + defaultBaselineMean : { # TODO - this should be the pedestal + full : @local::STM.HPGe.defaultBaselineMean + suppressed : 0 + } + defaultBaselineSD : { # TODO - this should be NoiseSD converted to ADC bins + full: @local::STM.HPGe.defaultBaselineSD + suppressed: 0 + } + STMMWDDigiTag : @local::DataProducts.PeakEnergies.LaBr } } } -# MakeTree_module.cc parameters -SimplifyStage1Data : { - StepPointMCsTag : "compressDetStepMCsSTM:" - StepPointMCsTag1809 : "compressDetStepMCsSTM:virtualdetector" - SimParticlemvTag : "compressDetStepMCsSTM:" - VirtualDetectorID : 101 +# ROOTAnalysisDump.fcl parameters +ROOTAnalysisDump : { + Stage1 : { + StepPointMCsTag : @local::DataProducts.Stage1.StepPointMCs + SimParticlemvTag : @local::DataProducts.Stage1.SimParticles + VirtualDetectorID : @local::DataProducts.Stage1.EndVirtualdetectorID + } + Stage2 : { + StepPointMCsTagVirtualdetector : @local::DataProducts.Stage2.StepPointMCs.Virtualdetector + StepPointMCsTagSTMDet : @local::DataProducts.Stage2.StepPointMCs.Detector + SimParticlemvTag : @local::DataProducts.Stage2.SimParticles + VirtualDetectorID : { + HPGe : @local::DataProducts.Stage2.EndVirtualdetectorID.HPGe + LaBr : @local::DataProducts.Stage2.EndVirtualdetectorID.LaBr + } + DetectorName : { + HPGe : "HPGe" + LaBr : "LaBr" + } + } + HPGeAnalysis: { + MicrospillHPGeWaveformTag : @local::DataProducts.Waveforms.HPGe.uSpill + ZSWaveformTag : @local::DataProducts.Waveforms.HPGe.ZS + MWDResultsTag : @local::DataProducts.PeakEnergies.HPGe + } + LaBrAnalysis: { + MicrospillHPGeWaveformTag : @local::DataProducts.Waveforms.LaBr.uSpill + ZSWaveformTag : @local::DataProducts.Waveforms.LaBr.ZS + MWDResultsTag : @local::DataProducts.PeakEnergies.LaBr + } + SimParticleVirtualDetectorBacktrace : { + StepPointMCsTag : { + Stage1 : @local::DataProducts.Stage1.StepPointMCs + Stage2 : @local::DataProducts.Stage2.StepPointMCs + } + StartVirtualDetectorID : { + Stage1 : @local::DataProducts.Stage1.EndVirtualdetectorID + Stage2 : { + HPGe : @local::DataProducts.Stage2.EndVirtualdetectorID.HPGe + LaBr : @local::DataProducts.Stage2.EndVirtualdetectorID.LaBr + } + } + TraceVirtualdetectorIds : { + Stage1: [101, 100, 10, 9, 8] + Stage2: [90, 101, 100, 10, 9, 8] + } + } } -SimplifyStage2Data : { - StepPointMCsTagSTMDet : "compressSTMDet:STMDet" - StepPointMCsTagVirtualdetector : "compressSTMDet:virtualdetector" - SimParticlemvTag : "compressSTMDet:" - VirtualDetectorID : { - HPGe : 90 - LaBr : 89 - } - DetectorName : { - HPGe : "HPGe" - LaBr : "LaBr" + +# For analysis studies, you should only need to edit lines between the "edit lines" comments +######################################## Start of edit lines ######################################## + +# Mixing parameters +MixSTMEvents : { + extendedMean2BB : 3.93e7 # Copied from Production/JobConfig/mixing/TwoBB.fcl + cutMax2BB : 2.36e8 # Copied from Production/JobConfig/mixing/TwoBB.fcl + extendedMean1BB : 1.58e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + cutMax1BB : 9.48e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + SDF : 0.6 # Copied from Production/JobConfig/mixing/OneBB.fcl + nPOTs : 26919873604557116 # Effective number of POTs, correcting for the resampling factor + nMicroSpills : 690253169 + meanEventsPerPOTFactors : { + EleBeamCat : 5.465876332164181e-11 + MuBeamCat : 7.887481312841522e-13 + TargetStopsCat1809 : 3.114418308224206e-16 } } @@ -174,25 +255,80 @@ STMDAQParameters : { LaBr : 370.370370370 # MSPS - copied from Offline/STMConditions/fcl/prolog.fcl } } + +# HPGe waveform generation parameters HPGeDigitization : { - EnergyPerADCBin : 0.5 # keV/bin + EnergyPerADCBin : 1 # keV/bin PreamplifierNoiseSD : 0.32 # mV DecayConstant : 50 # us - resetEventNumber : 31858 - microspillBufferLengthCount : 1000 + resetEventNumber : @local::DataProducts.nMicrospillsPerSpill + microspillBufferLengthCount : 2 concatenation : { - STMWaveformDigisTag : "DigiHPGe:" - nMerge : 31858 - filterTag : "concatenateWaveformsHPGe:" + STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill + nMerge : @local::DataProducts.nMicrospillsPerSpill + filterTag : @local::DataProducts.Waveforms.HPGe.concatenated } } -# Normalization parameter determination studies -# Detector efficiency simulation +# LaBr waveform generation parameters - not implemented yet +LaBrDigitization : { + EnergyPerADCBin : @nil # keV/bin + PreamplifierNoiseSD : @nil # mV + DecayConstant : @nil # us + resetEventNumber : @nil + microspillBufferLengthCount : @nil + concatenation : { + STMWaveformDigisTag : @local::DataProducts.Waveforms.LaBr.uSpill + nMerge : @nil + filterTag : @local::DataProducts.Waveforms.LaBr.concatenated + } +} + +# Detector efficiency simulation parameters Efficiency : { NPhotons : 1e7 - PhotonEnergy : 1.809 # MeV - OutputFilename : "Efficiency.root" + PhotonEnergy : @nil # MeV + OutputFilename : "nts.owner.Absorber.Spectrum.sequencer.root" + ST : { + x : -3904.0 # mm + y : 0.0 # mm + z : 627.0 # mm + } + FromSTToSSCAperture: { + HPGe: { + x : -40.6 # mm + y : 0.0 # mm + z : 39777 # mm + } + LaBr: { + x : 40.6 # mm + y : 0.0 # mm + z : 39777 # mm + } + } + FromSTToDet: { # ST at (-3904, 0, 627) + HPGe: { + x : -40.6 # mm + y : 0.0 # mm + z : 40072.1 # mm + } + LaBr: { + x : 40.6 # mm + y : 0.0 # mm + z : 40264.5 # mm + } + } +} +DigitizationTester: { # For firing photons into the detector to validate the behaviour of the relevant algorithms + HPGe: { + x : -3958.55 # mm + y : 0.0 # mm + z : 40584.95 # mm + deltax : -1.0 # MeV/c + deltay : 0.0 # MeV/c + deltaz : 1.0 # MeV/c + } } +######################################### End of edit lines ######################################### END_PROLOG diff --git a/STMMC/src/ConcatenationFilter_module.cc b/STMMC/src/ConcatenationFilter_module.cc index ed360dee97..10d6712377 100644 --- a/STMMC/src/ConcatenationFilter_module.cc +++ b/STMMC/src/ConcatenationFilter_module.cc @@ -58,6 +58,7 @@ namespace mu2e { void ConcatenationFilter::endJob() { mf::LogInfo log("ConcatenationFilter summary"); + log << "\n"; log << "==========ConcatenationFilter summary==========\n"; log << std::left << std::setw(25) << "\tNo. input events:" << inputEvents << "\n"; log << std::left << std::setw(25) << "\tNo. kept events:" << keptEvents << "\n"; diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index f084862a4e..19b2584492 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -1,6 +1,11 @@ // Simulates the electronics response of the HPGe detector. Simulates the pulse height, decay tail, and ADC digitization. Generates one STMWaveformDigi per micropulse. // Model based heavily on example provided in docDb 43617 // See docDb 51487 for full documentation +// Remaining TODOs for future development: +// - Include the commented out includes correctly +// - Implement the makeADCPlot functionality correctly +// - Implement different sampling frequencies correctly +// - Remove the implementation of stepPositionTolerance, it is a hack to address a position resolution I did not have time to fully address // Original author: Pawel Plesniak // stdlib includes @@ -45,6 +50,7 @@ // ROOT includes #include "art_root_io/TFileService.h" #include "TTree.h" +#include "TGraph.h" namespace mu2e { @@ -62,6 +68,7 @@ namespace mu2e { fhicl::Atom risingEdgeDecayConstant{ Name("risingEdgeDecayConstant"), Comment("Rising edge decay time [us]")}; fhicl::OptionalAtom microspillBufferLengthCount{ Name("microspillBufferLengthCount"), Comment("Number of microspills to buffer ahead for, in number of microspills")}; fhicl::OptionalAtom makeTTree{ Name("makeTTree"), Comment("Controls whether to make the TTree with branches chargeCollected, chargeDecayed, ADC, eventId, time")}; + fhicl::OptionalAtom makeADCPlot{ Name("makeADCPlot"), Comment("Controls whether to make a plot of the generated ADC values")}; fhicl::OptionalAtom timeOffset{ Name("timeOffset"), Comment("For debugging, adds the named time offset in [ns], used for testing analysis algorithms")}; fhicl::OptionalAtom resetEventNumber{ Name("resetEventNumber"), Comment("Simulates the off-spill period by resetting the inter-event last decayed charge to zero")}; fhicl::OptionalAtom verbosityLevel{ Name("verbosityLevel"), Comment("Controls verbosity")}; @@ -75,6 +82,7 @@ namespace mu2e { void decayCharge(); void addNoise(); void digitize(); + void endJob(); // fhicl variables art::ProductToken StepPointMCsTokenEle, StepPointMCsTokenMu, StepPointMCsToken1809; // Token of StepPointMCs in STMDet @@ -83,6 +91,7 @@ namespace mu2e { double noiseSD = 0; // Standard deviation of ADC noise [mV] double risingEdgeDecayConstant = 0; // [us] bool makeTTree = false; // Controls whether an analysis TTree is made + bool makeADCPlot = false; // Controls whether the ADC graph is made double timeOffset = 0.0; // Used for debugging [ns] uint resetEventNumber = 0; // Event ID at which to reset the charge amplitude int verbosityLevel = 0; // How much output to generate @@ -91,6 +100,7 @@ namespace mu2e { const double feedbackCapacitance = 1e-12; // [Farads] const double epsilonGe = 2.96; // Energy required to generate an eh pair in Ge at 77K [eV] const double micropulseTime = 1695.0; // [ns] + const double preamplifier_rise_time = 80.0; // [ns] // Define physics constants const double _e = 1.602176634e-19; // Electric charge constant [Coulombs] @@ -99,7 +109,8 @@ namespace mu2e { // ADC variables double chargeToADC = 0; // Conversion factor from charge built in capacitor to ADC determined voltage, multiply by this value to get from charge built to ADC voltage. - uint nADCs = 0; // Number of ADC values in an event + uint nADCs = 0; // Number of ADC values in an event, can change to account for timing variations + uint nADCs_init = 0; // Number of ADC values in an event, is fixed const int16_t ADCMax = static_cast((-1 * std::pow(2, 15)) + 1); // Maximum ADC value, power is 15 not 16 as using int16_t not uint16_t double ADC = 0; // iterator variable uint32_t eventTimeBuffer = 0; // Multiple of event ids to store @@ -137,7 +148,8 @@ namespace mu2e { double electronTravelTime = 0, holeTravelTime = 0; // Drift times [ns] uint32_t electronTravelTimeSteps = 0, holeTravelTimeSteps = 0; // Drift times [steps] double decayExp = 0; // Amount of decay with each tADC - double lastEventEndDecayedCharge = 0; // Carry over for starting new microspill waveforms [charge carrier pairs] + double lastEventEndDecayedCharge = 0; // Carry over the charge stored in the crystal from the previous event + double lastEventEndIntegratedCharge = 0; // Carry over for integrated charge stored in the preamplifier int microspillBufferLengthCount = 0; // Buffer to store the charge deposits that are allocated to this event but happen after the microspill ends e.g. 844keV const int defaultMicrospillBufferLengthCount = 2; // Default value for the microspill buffer length @@ -155,6 +167,10 @@ namespace mu2e { std::vector _chargeCarryOver; // Temporary buffer that will store _chargeCollected over the course of the next event std::vector _adcs; // Buffer for storing the ADC values to put into the STMWaveformDigi + // Debugging variables + double waveform_time = 0.0; // Time of start of waveform for debugging. TODO - REMOVE ME! + int32_t kept_events = 0, dropped_events = 0; // Counters for debugging + // Offline utilities // TODO: include the prodition to get the sampling frequency // mu2e::STMChannel::enum_type _HPGeChannel = static_cast(1); @@ -172,19 +188,19 @@ namespace mu2e { noiseSD(conf().noiseSD()), risingEdgeDecayConstant(conf().risingEdgeDecayConstant()) { produces(); - if (defaultMicrospillBufferLengthCount < 2) - throw cet::exception("RANGE", "defaultMicrospillBufferLengthCount has to be more than 1\n"); - crystalCentrePosition.set(crystalCentreX, crystalCentreY, crystalCentreZ); tADC = 1e3/fADC; // 1e3 converts [us] to [ns] as fADC is in [MHz] // Assign optional variables microspillBufferLengthCount = conf().microspillBufferLengthCount() ? *(conf().microspillBufferLengthCount()) : defaultMicrospillBufferLengthCount; + if (defaultMicrospillBufferLengthCount < 1) + throw cet::exception("RANGE", "defaultMicrospillBufferLengthCount has to be at least 1!\n"); verbosityLevel = conf().verbosityLevel() ? *(conf().verbosityLevel()) : 0; // Determine the number of ADC values in each STMWaveformDigi. Increase the number by one due to truncation. At 320MHz, this will be 543 ADC values per microbunch double _nADCs = (micropulseTime/tADC) + 1; - nADCs = (int) _nADCs; + nADCs = (uint) _nADCs; + nADCs_init = nADCs; _charge.insert(_charge.begin(), nADCs * microspillBufferLengthCount, 0.); _chargeCollected.insert(_chargeCollected.begin(), nADCs * microspillBufferLengthCount, 0.); _chargeCarryOver.insert(_chargeCarryOver.begin(), nADCs * (microspillBufferLengthCount - 1), 0.); @@ -205,7 +221,7 @@ namespace mu2e { // Assign the approrpiate time offset timeOffset = conf().timeOffset() ? *(conf().timeOffset()) : 0.0; - // Assign TTrees + // Assign optional ROOT parameters makeTTree = conf().makeTTree() ? *(conf().makeTTree()) : false; if (makeTTree) { art::ServiceHandle tfs; @@ -216,13 +232,22 @@ namespace mu2e { ttree->Branch("eventId", &eventId, "eventId/i"); ttree->Branch("time", &time, "time/i"); }; + makeADCPlot = conf().makeADCPlot() ? *(conf().makeADCPlot()) : false; + if (makeADCPlot) { + throw cet::exception("IMPLEMENT", "makeADCPlot functionality not yet implemented correctly\n"); + art::ServiceHandle tfs; + TTree* adcTree = tfs->make("adcTree", "HPGeWaveformsFromStepPointMCs ADC Tree"); + adcTree->Branch("ADC", &ADC, "ADC/D"); + }; + // Assign optional variables resetEventNumber = conf().resetEventNumber() ? *(conf().resetEventNumber()) : 0; }; void HPGeWaveformsFromStepPointMCs::beginJob() { if (verbosityLevel) { - std::cout << "STM HPGe digitization parameters" << std::endl; + std::cout << std::endl; + std::cout << "=======================================STM HPGe digitization parameters=======================================" << std::endl; std::cout << "\tInput parameters" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "fAD [MHz]" << fADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "EnergyPerADCBin [keV/bin]" << ADCToEnergy << std::endl; @@ -230,12 +255,14 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "risingEdgeDecayConstant [us]" << risingEdgeDecayConstant << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "microspillBufferLengthCount" << microspillBufferLengthCount << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "makeTTree" << makeTTree << std::endl; + std::cout << std::left << "\t\t" << std::setw(60) << "makeADCPlot" << makeADCPlot << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "timeOffset [ns]" << timeOffset << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "resetEventNumber" << resetEventNumber << std::endl; std::cout << "\tDerived parameters: " << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "tADC [ns]" << tADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "nADCs" << nADCs << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "NoiseSD [charge carriers]" << noiseSD << std::endl; + std::cout << std::left << "\t\t" << std::setw(60) << "NoiseSD [ADC bins]" << noiseSD * chargeToADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "chargeToADC [bin/charge carrier]" << chargeToADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range [V]" << "[+1, -1]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [V]" << "[0, -1]" << std::endl; @@ -243,42 +270,80 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [charge carriers]" << "[0, " << ADCMax/chargeToADC << "]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [C]" << "[0, " << ADCMax * _e/chargeToADC << "]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Energy range [keV]" << "[0, " << -1*ADCMax * ADCToEnergy << "]" << std::endl; + std::cout << "==============================================================================================================" << std::endl; std::cout << std::endl; // buffer line }; }; void HPGeWaveformsFromStepPointMCs::produce(art::Event& event) { eventId = event.id().event(); + // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO + // Create the STMWaveformDigi and insert all the relevant attributes + // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times + if (std::abs(fADC - 320) > 1e-12) + throw cet::exception("RANGE") << "Currently only fADC of 320 MHz is supported (got " << fADC << " MHz)\n"; + nADCs = nADCs_init; + + // Assign the variables for time and number of ADCs + eventTimeBuffer = eventId % 5; + if (!(eventTimeBuffer == 0 || eventTimeBuffer == 3)) + nADCs++; + + // Update the last event decayed charge before the noise so the noise isn't added twice + if (resetEventNumber != 0 && eventId == resetEventNumber) { + lastEventEndDecayedCharge = 0.0; + lastEventEndIntegratedCharge = 0.0; + std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); + } + + // Update the parameters to carry over to the next event + _chargeCollected.clear(); + _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); + _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); + _chargeCarryOver.clear(); + _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); + + // Clear previous buffer vectors + _chargeDecayed.clear(); + _chargeDecayed.insert(_chargeDecayed.end(), nADCs, 0); + _adcs.clear(); + _adcs.insert(_adcs.end(), nADCs, 0); + + // Reset waveform time for debugging + waveform_time = 0.0; + // Get the hits in the detector std::vector StepsEle = event.getProduct(StepPointMCsTokenEle); std::vector StepsMu = event.getProduct(StepPointMCsTokenMu); std::vector Steps1809 = event.getProduct(StepPointMCsToken1809); - // Add a collection of charge depositions to _charge - for(const StepPointMC& step : StepsEle){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; - for(const StepPointMC& step : StepsMu){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; - for(const StepPointMC& step : Steps1809){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; + // If there are no hits, produce an empty STMWaveformDigiCollection and return + bool hasHits = (!StepsEle.empty() || !StepsMu.empty() || !Steps1809.empty()); + + if (hasHits) { + kept_events++; + // Add a collection of charge depositions to _charge + for(const StepPointMC& step : StepsEle) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + for(const StepPointMC& step : StepsMu) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + for(const StepPointMC& step : Steps1809) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + } else { + dropped_events++; + // We do nothing here. _chargeCollected is already 0s from initialization. + // We just proceed to decayCharge(). + } // Decay all of the collected charges decayCharge(); - // Update the last event decayed charge before the noise so the noise isn't added twice - if (resetEventNumber != 0 && eventId == resetEventNumber) { - lastEventEndDecayedCharge = 0; - std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); - } - else - lastEventEndDecayedCharge = _chargeDecayed.back(); - // Add preamplifier electronics noise with SD defined in noiseSD addNoise(); @@ -291,21 +356,13 @@ namespace mu2e { throw cet::exception("LogicError", "ADC values too high!"); }; - // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO - // Create the STMWaveformDigi and insert all the relevant attributes - // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times - eventTimeBuffer = eventId % 5; - if (eventTimeBuffer == 0 || (eventId % 3) == 0) - eventTime += nADCs + 1; - else - eventTime += nADCs; STMWaveformDigi _waveformDigi(eventTime, _adcs); std::unique_ptr outputDigis(new STMWaveformDigiCollection); outputDigis->emplace_back(_waveformDigi); // Make the ttree if appropriate if (makeTTree) { - time = _waveformDigi.trigTimeOffset() - 1; + time = eventTime; for (uint i = 0; i < nADCs; i++) { chargeCollected = _chargeCollected[i]; chargeDecayed = _chargeDecayed[i]; @@ -315,18 +372,27 @@ namespace mu2e { }; }; - // Update the parameters to carry over to the next event - _chargeCarryOver.clear(); - _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); - _chargeCollected.clear(); - _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); - _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); - // Clear previous buffer vectors - std::fill(_chargeDecayed.begin(), _chargeDecayed.end(), 0); - std::fill(_adcs.begin(), _adcs.end(), 0); + // Make the ADC plot if appropriate + if (makeADCPlot) { + // NOTE - THIS DOES NOT CURRENTLY FUNCTION, AND I DO NOT HAVE TIME TO FIX IT RIGHT NOW + std::stringstream histsuffix; + histsuffix.str(""); + histsuffix << "_evt" << event.event(); + + art::ServiceHandle tfs; + std::vector _adc_double = std::vector(_adcs.begin(), _adcs.end()); + TGraph* g_adcs = tfs->make(nADCs, _adc_double.data()); + g_adcs->SetName(("g_adcs"+histsuffix.str()).c_str()); + // g_adcs->SetTitle("ADC Values"); + // for (uint i = 0; i < nADCs; i++) + // g_adcs->SetPoint(i, i, _adcs[i]); + }; // Add the STMWaveformDigi to the event event.put(std::move(outputDigis)); + + // Update the event time for the next waveform + eventTime += nADCs; return; }; @@ -436,7 +502,8 @@ namespace mu2e { }; // Allocate the rest of the charge for one more entry for the continuity - _charge[tIndex] = N_ehPairs; + for (uint i = tIndex; i < _charge.size(); i++) + _charge[i] = _charge[tIndex - 1]; tIndex++; // Update _chargeCollected. First case is treated separately as there is no charge deposited in the previous step @@ -453,10 +520,28 @@ namespace mu2e { }; void HPGeWaveformsFromStepPointMCs::decayCharge() { - _chargeDecayed[0] = lastEventEndDecayedCharge * decayExp + _chargeCollected[0]; - for (uint t = 1; t < nADCs; t++) - _chargeDecayed[t] = _chargeDecayed[t-1] * decayExp + _chargeCollected[t]; - return; + // Convert time constants to samples based on 320MHz clock, 80ns is the approximate ORTEC GMX preamplifier rise time + const double tau_rise_samples = preamplifier_rise_time / tADC; // Convert this to ADC steps + const double riseExp = 1.0 / (1.0 + tau_rise_samples); + + // Carry over the current in the last time step + double currentDecayed = lastEventEndDecayedCharge; + + // Holiding variable + double currentIntegrated = lastEventEndIntegratedCharge; + + for (uint t = 0; t < nADCs; t++) { + // Discharge the capacitor, add the next contribution of the step + currentDecayed = currentDecayed * decayExp + _chargeCollected[t]; + + // Smooth out the waveform accounting for the finite bandwidth of the preamplifier + currentIntegrated = riseExp * currentDecayed + (1.0 - riseExp) * currentIntegrated; + + // Store the results + _chargeDecayed[t] = currentIntegrated; + } + lastEventEndDecayedCharge = currentDecayed; + lastEventEndIntegratedCharge = currentIntegrated; }; void HPGeWaveformsFromStepPointMCs::addNoise() { @@ -476,10 +561,19 @@ namespace mu2e { // Convert the charge deposition to ADC voltage output. for (uint i = 0; i < nADCs; i++) { ADC = _chargeDecayed[i] * chargeToADC; - _adcs[i] = ADC > ADCMax ? static_cast(std::round(ADC)) : ADCMax; + _adcs[i] = (ADC > ADCMax) ? static_cast(std::round(ADC)) : ADCMax; }; return; }; + + void HPGeWaveformsFromStepPointMCs::endJob() { + mf::LogInfo log("HPGeWaevfrmsFromStepPointMCs"); + log << "\n"; + log << "=====HPGeWaveformsFromStepPointMCs summary=====\n"; + log << std::left << std::setw(25) << "\tNo. kept events: " << kept_events << "\n"; + log << std::left << std::setw(25) << "\tNo. discarded events: " << dropped_events << "\n"; + log << "===============================================\n"; + }; }; // namespace mu2e DEFINE_ART_MODULE(mu2e::HPGeWaveformsFromStepPointMCs) diff --git a/STMMC/src/STMResamplingProducer_module.cc b/STMMC/src/STMResamplingProducer_module.cc index 02c83d003f..30137f61df 100644 --- a/STMMC/src/STMResamplingProducer_module.cc +++ b/STMMC/src/STMResamplingProducer_module.cc @@ -44,7 +44,7 @@ namespace mu2e { art::EDProducer{conf}, StepPointMCsToken(consumes(conf().stepPointMCsTag())), virtualDetectorID(conf().virtualDetectorID()) { - produces(); + produces("virtualdetector"); }; void STMResamplingProducer::produce(art::Event& event) { @@ -63,7 +63,7 @@ namespace mu2e { // Update counter includedStepPointMCs += outputStepPointMCs->size(); // Add the new data products to the event - event.put(std::move(outputStepPointMCs)); + event.put(std::move(outputStepPointMCs), "virtualdetector"); return; }; diff --git a/STMMC/src/SimParticleAndVDBacktrace_module.cc b/STMMC/src/SimParticleAndVDBacktrace_module.cc new file mode 100644 index 0000000000..82f4538588 --- /dev/null +++ b/STMMC/src/SimParticleAndVDBacktrace_module.cc @@ -0,0 +1,263 @@ +// Adapted from ReadVirtualDetector_module.cc +// Generates a large ROOT TTree with SimParticle and StepPointMC information. For each SimParticle that has a StepPointMC in the specified virtual detector, gets its +// - parentCounter - 0 for the original particle, 1 for its parent, 2 for its grandparent, etc. This is traced all the way back to the POT +// - particleId - SimParticle ID (from GEANT) +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// For each chosen virtual detector ID to trace back to, also gets the following StepPointMC information (0 if the SimParticle did not have a StepPointMC in that virtual detector): +// - VD__x [mm] - x position of the StepPointMC in that virtual detector +// - VD__y [mm] - y position of the StepPointMC in that virtual detector +// - VD__z [mm] - z position of the StepPointMC in that virtual detector +// - VD__px [MeV/c] - x momentum of the StepPointMC in that virtual detector +// - VD__py [MeV/c] - y momentum of the StepPointMC in that virtual detector +// - VD__pz [MeV/c] - z momentum of the StepPointMC in that virtual detector +// - VD__time [ns] - time of the StepPointMC in that virtual detector +// - VD__Ekin [MeV] - kinetic energy of the StepPointMC in that virtual detector +// - VD__pdgId - PDG ID of the SimParticle at that StepPointMC +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleAndVDBacktrace : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom StartVirtualDetectorId{Name("StartVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::Sequence TraceVirtualdetectorIds{Name("TraceVirtualdetectorIds"), Comment("IDs of the virtual detectors to trace back to the origin")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleAndVDBacktrace(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const art::Ptr &particle, const std::vector steps); + void addParentToTree(const art::Ptr &particle, const std::vector steps); + void traceSteps(const art::Ptr &particle, const std::vector steps); + void clear_virtualdetector_vectors(); + void populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index); + private: + art::ProductToken StepPointMCsToken; + unsigned StartVirtualDetectorId = 0; + std::vector TraceVirtualdetectorIds; + int consecutiveEmptyFileThreshold = 0; + + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, nTraceVirtualdetectorIds = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + std::vector VD_x, VD_y, VD_z, VD_px, VD_py, VD_pz, VD_time, VD_Ekin; + std::vector VD_pdgId; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + std::vector*> VD_vectors; + }; + + SimParticleAndVDBacktrace::SimParticleAndVDBacktrace(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + StartVirtualDetectorId(conf().StartVirtualDetectorId()), + TraceVirtualdetectorIds(conf().TraceVirtualdetectorIds()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + if (std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), StartVirtualDetectorId) == TraceVirtualdetectorIds.end()) { + TraceVirtualdetectorIds.push_back(StartVirtualDetectorId); + }; + nTraceVirtualdetectorIds = TraceVirtualdetectorIds.size(); + VD_vectors = {&VD_x, &VD_y, &VD_z, &VD_px, &VD_py, &VD_pz, &VD_time, &VD_Ekin}; + for (auto vec : VD_vectors) { + vec->resize(nTraceVirtualdetectorIds, 0.0); + }; + VD_pdgId.resize(nTraceVirtualdetectorIds, 0); + + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + for (size_t i=0; iBranch(Form("VD_%u_x", TraceVirtualdetectorIds[i]), &VD_x[i], Form("VD_%u_x/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_y", TraceVirtualdetectorIds[i]), &VD_y[i], Form("VD_%u_y/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_z", TraceVirtualdetectorIds[i]), &VD_z[i], Form("VD_%u_z/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_px", TraceVirtualdetectorIds[i]), &VD_px[i], Form("VD_%u_px/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_py", TraceVirtualdetectorIds[i]), &VD_py[i], Form("VD_%u_py/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_pz", TraceVirtualdetectorIds[i]), &VD_pz[i], Form("VD_%u_pz/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_time", TraceVirtualdetectorIds[i]), &VD_time[i], Form("VD_%u_time/D", TraceVirtualdetectorIds[i])); // ns + ttree->Branch(Form("VD_%u_Ekin", TraceVirtualdetectorIds[i]), &VD_Ekin[i], Form("VD_%u_Ekin/D", TraceVirtualdetectorIds[i])); // MeV + ttree->Branch(Form("VD_%u_pdgId", TraceVirtualdetectorIds[i]), &VD_pdgId[i], Form("VD_%u_pdgId/I", TraceVirtualdetectorIds[i])); + }; + }; + + void SimParticleAndVDBacktrace::populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index) { + const CLHEP::Hep3Vector momentum = step_it.momentum(); + const CLHEP::Hep3Vector position = step_it.position(); + const double mass = pdt->particle(step_it.simParticle()->pdgId()).mass(); + const double EKin = std::sqrt(momentum.mag2() + mass * mass) - mass; + VD_x[index] = position.x(); + VD_y[index] = position.y(); + VD_z[index] = position.z(); + VD_px[index] = momentum.x(); + VD_py[index] = momentum.y(); + VD_pz[index] = momentum.z(); + VD_time[index] = step_it.time(); + VD_Ekin[index] = EKin; + VD_pdgId[index] = step_it.simParticle()->pdgId(); + }; + + void SimParticleAndVDBacktrace::traceSteps(const art::Ptr &particle, const std::vector steps) { + unsigned step_virtualdetectorId = 0, index = 0; + for (const StepPointMC & step_it : steps) { + if (step_it.simParticle()->id() != particle->id()) continue; + step_virtualdetectorId = step_it.virtualDetectorId(); + auto it = std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), step_virtualdetectorId); + if (it == TraceVirtualdetectorIds.end()) continue; + index = std::distance(TraceVirtualdetectorIds.begin(), it); + populate_virtualdetector_vector(step_it, index); + }; + }; + + void SimParticleAndVDBacktrace::clear_virtualdetector_vectors() { + for (auto & vec : VD_vectors) { + std::fill(vec->begin(), vec->end(), 0.0); + }; + std::fill(VD_pdgId.begin(), VD_pdgId.end(), 0); + return; + }; + + void SimParticleAndVDBacktrace::addToTree(const art::Ptr &particle, const std::vector steps) { + // Assigne all the SimParticle variables to the tree variables + parentCounter = 0; + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + + // Now trace the parents + if (!particle->hasParent()) return; + + // Require separate definitions of addParentToTree because of data accessors + parent = particle->parent(); + addParentToTree(parent, steps); + + while (parent->hasParent()) { + parent = parent->parent(); + addParentToTree(parent, steps); + }; + }; + + void SimParticleAndVDBacktrace::addParentToTree(const art::Ptr &particle, const std::vector steps) { + // Increment the parent counter + parentCounter++; + + // Assign all the SimParticle variables to the tree variables + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + }; + + void SimParticleAndVDBacktrace::analyze(const art::Event& event) { + // Get the data products + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, are you sure you have the correct data product name?\n"); + }; + return; + }; + const std::vector StepPointMCs = *stepHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != StartVirtualDetectorId) continue; + + // Get the associated particle + const art::Ptr particle = step.simParticle(); + addToTree(particle, StepPointMCs); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleAndVDBacktrace) diff --git a/STMMC/src/SimParticleDump_module.cc b/STMMC/src/SimParticleDump_module.cc new file mode 100644 index 0000000000..48cf595cf6 --- /dev/null +++ b/STMMC/src/SimParticleDump_module.cc @@ -0,0 +1,182 @@ +// Adapted from ReadVirtualDetector_module.cc +// For StepPointMCs in virtualdetectors, generates a TTree showing the provenance of the generated SimParticles with the following branches +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// TODO before upload - validate that the code works and generates the same particle trace as the SimParticleAndVDBacktrace module + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleDump : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom SimParticlemvTag{Name("SimParticlemvTag"), Comment("Tag identifying the SimParticlemv")}; + fhicl::Atom FilterVirtualDetectorId{Name("FilterVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleDump(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const SimParticle& particle, TTree* ttree); + void addParentToTree(art::Ptr &_particle, TTree* ttree); + private: + art::ProductToken StepPointMCsToken; + art::ProductToken SimParticlemvToken; + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + unsigned virtualdetectorId = 0, FilterVirtualDetectorId = 0; + bool forward = false; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + }; + + SimParticleDump::SimParticleDump(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + SimParticlemvToken(consumes(conf().SimParticlemvTag())), + FilterVirtualDetectorId(conf().FilterVirtualDetectorId()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + }; + + void SimParticleDump::addToTree(const SimParticle& particle, TTree* ttree) { + parentCounter = 0; + particleId = particle.id(); + pdgId = particle.pdgId(); + creationCode = particle.creationCode(); + startPosition = particle.startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle.startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle.startGlobalTime(); + Ekin = particle.endKineticEnergy(); + ttree->Fill(); + + if (!particle.hasParent()) { + return; + }; + parentCounter++; + parent = particle.parent(); + addParentToTree(parent, ttree); + while (parent->hasParent()) { + parentCounter++; + parent = parent->parent(); + addParentToTree(parent, ttree); + }; + }; + + void SimParticleDump::addParentToTree(art::Ptr &_particle, TTree* ttree) { + particleId = _particle->id(); + pdgId = _particle->pdgId(); + creationCode = _particle->creationCode(); + startPosition = _particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = _particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = _particle->startGlobalTime(); + Ekin = _particle->endKineticEnergy(); + ttree->Fill(); + }; + + void SimParticleDump::analyze(const art::Event& event) { + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + auto simHandle = event.getHandle< SimParticleCollection >(SimParticlemvToken); + if (!simHandle || simHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, stopping the job"); + }; + + auto const& StepPointMCs = *stepHandle; + auto const& SimParticles = *simHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != FilterVirtualDetectorId) continue; + + // Get the associated particle + const SimParticle& particle = SimParticles.at(step.trackId()); + addToTree(particle, ttree); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleDump) diff --git a/STMMC/src/VirtualDetectorTree_module.cc b/STMMC/src/VirtualDetectorTree_module.cc index 6156d1fc6b..a976156ce1 100644 --- a/STMMC/src/VirtualDetectorTree_module.cc +++ b/STMMC/src/VirtualDetectorTree_module.cc @@ -19,6 +19,7 @@ // fhicl includes #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" // message handling #include "messagefacility/MessageLogger/MessageLogger.h" @@ -45,6 +46,7 @@ namespace mu2e { struct Config { fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; fhicl::Atom SimParticlemvTag{Name("SimParticlemvTag"), Comment("Tag identifying the SimParticlemv")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; }; using Parameters = art::EDAnalyzer::Table; explicit VirtualDetectorTree(const Parameters& conf); @@ -54,17 +56,19 @@ namespace mu2e { art::ProductToken StepPointMCsToken; art::ProductToken SimParticlemvToken; GlobalConstantsHandle pdt; - int pdgId = 0; - double x = 0.0, y = 0.0, z = 0.0, mass = 0.0, E = 0.0, time = 0.0; + int pdgId = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0, simParticleId=0; + double x = 0.0, y = 0.0, z = 0.0, mass = 0.0, Ekin = 0.0, Etot = 0.0, time = 0.0, p = 0.0, p2 = 0.0, px = 0.0, py = 0.0, pz = 0.0; VolumeId_type virtualdetectorId = 0; TTree* ttree; std::map pdgIds; // + uint simParticleIdKey = 0; }; VirtualDetectorTree::VirtualDetectorTree(const Parameters& conf) : art::EDAnalyzer(conf), StepPointMCsToken(consumes(conf().StepPointMCsTag())), SimParticlemvToken(consumes(conf().SimParticlemvTag())) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; art::ServiceHandle tfs; ttree = tfs->make( "ttree", "Virtual Detectors ttree"); ttree->Branch("time", &time, "time/D"); // ns @@ -73,18 +77,36 @@ namespace mu2e { ttree->Branch("x", &x, "x/D"); // mm ttree->Branch("y", &y, "y/D"); // mm ttree->Branch("z", &z, "z/D"); // mm - ttree->Branch("E", &E, "E/D"); // MeV + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("p", &p, "p/D"); // mm + ttree->Branch("p2", &p2, "p2/D"); // mm + ttree->Branch("mass", &mass, "mass/D"); // MeV/c^2 + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + ttree->Branch("Etot", &Etot, "Etot/D"); // MeV + ttree->Branch("SimParticleId", &simParticleId, "SimParticleId/i"); }; void VirtualDetectorTree::analyze(const art::Event& event) { - // Get the data products from the event - auto const& StepPointMCs = event.getProduct(StepPointMCsToken); - if (StepPointMCs.empty()) - return; - auto const& SimParticles = event.getProduct(SimParticlemvToken); - if (SimParticles.empty()) + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; return; + } + // Try to get SimParticles (assuming it's a map or vector) + auto simHandle = event.getHandle< SimParticleCollection >(SimParticlemvToken); + if (!simHandle || simHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + } + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, stopping the job"); + } + auto const& StepPointMCs = *stepHandle; + auto const& SimParticles = *simHandle; + consecutiveEmptyFileCounter = 0; // Loop over all VD hits for (const StepPointMC& step : StepPointMCs) { // Get the associated particle @@ -97,9 +119,16 @@ namespace mu2e { x = step.position().x(); y = step.position().y(); z = step.position().z(); + px = step.momentum().x(); + py = step.momentum().y(); + pz = step.momentum().z(); + p2 = step.momentum().mag2(); + p = std::sqrt(p2); mass = pdt->particle(pdgId).mass(); - E = std::sqrt(step.momentum().mag2()+mass*mass)-mass; // Subtract the rest mass - if (E < 0) + Etot = std::sqrt(p2 + mass * mass); // Total energy + Ekin = Etot - mass; // Subtract the rest mass + simParticleId = particle.id().asInt(); + if (Ekin < 0) throw cet::exception("LogicError", "Energy is negative"); ttree->Fill(); diff --git a/STMReco/fcl/plotSTMWaveformDigis.fcl b/STMReco/fcl/plotSTMWaveformDigis.fcl index e6b1bca596..cd3e32eb7c 100644 --- a/STMReco/fcl/plotSTMWaveformDigis.fcl +++ b/STMReco/fcl/plotSTMWaveformDigis.fcl @@ -21,7 +21,7 @@ physics: { producers : { } filters : { } analyzers : { - + plotRawWaveformDigisHPGe : { module_type : PlotSTMWaveformDigis stmWaveformDigisTag : "makeSTMDigis:rawHPGe" @@ -29,7 +29,7 @@ physics: { verbosityLevel : 0 xAxis : "sample_number" } - + plotZSWaveformDigisHPGe : { module_type : PlotSTMWaveformDigis stmWaveformDigisTag : "makeSTMDigis:zsHPGe" @@ -53,7 +53,7 @@ physics: { verbosityLevel : 0 xAxis : "sample_number" } - + plotZSWaveformDigisLaBr : { module_type : PlotSTMWaveformDigis stmWaveformDigisTag : "makeSTMDigis:zsLaBr" @@ -69,7 +69,7 @@ physics: { verbosityLevel : 0 xAxis : "sample_number" } - + } # setup paths trigger_paths: [ ] diff --git a/STMReco/src/MakeSTMHits_module.cc b/STMReco/src/MakeSTMHits_module.cc index abde7f2e4b..a2c9fd27b6 100644 --- a/STMReco/src/MakeSTMHits_module.cc +++ b/STMReco/src/MakeSTMHits_module.cc @@ -59,7 +59,7 @@ namespace mu2e { produces(); } //Originally had _channel(STMChannel::LaBr) - + void MakeSTMHits::produce(art::Event& event) { // create output unique_ptr outputSTMHits(new STMHitCollection); diff --git a/STMReco/src/PlotSTMPHSpectrum_module.cc b/STMReco/src/PlotSTMPHSpectrum_module.cc index 4995dc4aef..e39790ccb2 100644 --- a/STMReco/src/PlotSTMPHSpectrum_module.cc +++ b/STMReco/src/PlotSTMPHSpectrum_module.cc @@ -46,7 +46,7 @@ namespace mu2e { TH2F* _twoDhist; //Histograms of Energy vs binned event int eventCount = 0; - + TH1D* _phSpectrum; art::ProductToken _stmPHDigisToken; }; @@ -61,15 +61,15 @@ namespace mu2e { // 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 + 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); diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 4bbe13112f..b032d52b4e 100644 --- a/STMReco/src/STMMovingWindowDeconvolution_module.cc +++ b/STMReco/src/STMMovingWindowDeconvolution_module.cc @@ -66,7 +66,7 @@ namespace mu2e { 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\"") - }; + }; }; using Parameters = art::EDProducer::Table; @@ -115,7 +115,15 @@ namespace mu2e { std::vector averaged_data; double sum = 0.0; // variable part of avarege() int count = 0; // counter variable - int16_t ph_energy = 0; + int16_t mwd_energy = 0; + bool _isFirstWaveform = true; + double _last_ADC = 0.0; + double _last_deconvolved = 0.0; + // Rolling history buffers to bridge the gap between chunks + std::deque _deconvolved_history; // Stores last M samples + std::deque _differentiated_history; // Stores last L samples + double _last_average = 0.0; // Stores the last average value from previous chunk + // Peak finding variables double baseline_mean = 0.0; @@ -140,7 +148,10 @@ namespace mu2e { uint eventId = 0, waveformID = 0; // Summary data - uint processedEvents = 0, processedWaveforms = 0, foundPeaks = 0; + uint processedEvents = 0, processedWaveforms = 0, foundPeaks = 0, processedADCs; + + // Debugging variables + double waveform_time_microspill_frame = 0.0; }; STMMovingWindowDeconvolution::STMMovingWindowDeconvolution(const Parameters& conf) : @@ -190,7 +201,8 @@ namespace mu2e { void STMMovingWindowDeconvolution::beginJob() { if (verbosityLevel) { - std::cout << "STM Moving Window Deconvolution" << std::endl; + std::cout << std::endl; + std::cout << "=======================================STM HPGe digitization parameters=======================================" << std::endl; std::cout << "\tAlgorithm parameters" << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "tau" << tau << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "M" << M << std::endl; @@ -200,6 +212,7 @@ namespace mu2e { std::cout << "\tChannel: " << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "Name" << channel.name() << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "ID" << static_cast(channel.id()) << std::endl; + std::cout << "==============================================================================================================" << std::endl; std::cout << std::endl; // buffer line }; }; @@ -222,9 +235,9 @@ namespace mu2e { //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 - STMPHDigi ph_digi(peak_times[i], -1 * ph_energy); // peak_heights are negative, make them positive here + mwd_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 * mwd_energy); // peak_heights are negative, make them positive here if (ph_digi.energy() < -100) throw cet::exception("logicError", "The peak height must be positive!"); @@ -395,11 +409,23 @@ namespace mu2e { void STMMovingWindowDeconvolution::endJob() { mf::LogInfo log("MWD summary"); + log << "\n"; log << "==================MWD summary==================\n"; log << std::left << std::setw(25) << "\tProcessed events: " << processedEvents << "\n"; log << std::left << std::setw(25) << "\tProcessed waveforms: " << processedWaveforms << "\n"; log << std::left << std::setw(25) << "\tFound peaks: " << foundPeaks << "\n"; + log << std::left << std::setw(25) << "\tnADCs: " << nADCs << "\n"; + log << std::left << std::setw(25) << "\tprocessedADCs: " << processedADCs << "\n"; log << "===============================================\n"; + if (verbosityLevel) { + log << "\n"; + log << "==================Parameters==================\n"; + log << std::left << std::setw(25) << "\tnsPerCt: " << nsPerCt << "\n"; + log << std::left << std::setw(25) << "\tPedestal: " << pedestal << "\n"; + log << std::left << std::setw(25) << "\tBaseline mean: " << baseline_mean << "\n"; + log << std::left << std::setw(25) << "\tBaseline std dev: " << baseline_stddev << "\n"; + log << "===============================================\n"; + } return; }; @@ -445,7 +471,7 @@ namespace mu2e { 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]); @@ -456,7 +482,7 @@ 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"); diff --git a/STMReco/test/plotDummyDigis.C b/STMReco/test/plotDummyDigis.C index 4e2accf708..facc6961f0 100644 --- a/STMReco/test/plotDummyDigis.C +++ b/STMReco/test/plotDummyDigis.C @@ -38,7 +38,7 @@ void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { //No need to redefine the histograms from before TCanvas* c2 = new TCanvas(); c2->SetTitle("STMWaveformDigis Super Imposed"); - + hRawWaveform->SetLineColor(kGreen); hRawWaveform->SetLineWidth(3); @@ -62,7 +62,7 @@ void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { //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()); @@ -74,4 +74,3 @@ void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { // c4->DrawClonePad(); } - diff --git a/TrkReco/src/TrkUtilities.cc b/TrkReco/src/TrkUtilities.cc index 27c4b025b2..034b6929c6 100644 --- a/TrkReco/src/TrkUtilities.cc +++ b/TrkReco/src/TrkUtilities.cc @@ -81,4 +81,3 @@ namespace mu2e { double gamma(double mass, double momentum) { return energy(mass,momentum)/mass; } } // TrkUtilities }// mu2e -