diff --git a/sbndcode/CMakeLists.txt b/sbndcode/CMakeLists.txt index 07f557b81..c13c94a16 100644 --- a/sbndcode/CMakeLists.txt +++ b/sbndcode/CMakeLists.txt @@ -47,3 +47,5 @@ add_subdirectory(ChannelMaps) add_subdirectory(SBNDCVN) add_subdirectory(TPCPMTBarycenterMatching) + +add_subdirectory(LightPropagationCorrection) \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/CMakeLists.txt b/sbndcode/LightPropagationCorrection/CMakeLists.txt new file mode 100644 index 000000000..a40b80774 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/CMakeLists.txt @@ -0,0 +1,51 @@ + +set ( + MODULE_LIBRARIES + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support + art_root_io::TFileService_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + + ROOT::Tree + ROOT::Core + + larsim::Utils + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + larcore::Geometry_Geometry_service + larcorealg::Geometry + lardataobj::RawData + lardataobj::RecoBase + lardataobj::MCBase + lardataobj::Simulation + nusimdata::SimulationBase + nug4::ParticleNavigation + + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_SBND + sbndaq_artdaq_core::sbndaq-artdaq-core_Obj_SBND + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + + sbndcode_OpDetReco_OpFlash_FlashFinder + + sbnobj::Common_Reco + sbnobj::SBND_CRT + sbnobj::SBND_Timing + sbnobj::SBND_OpFlashTiming + sbndcode_OpDetSim +) + +cet_build_plugin(LightPropagationCorrection art::module SOURCE LightPropagationCorrection_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(LightPropagationCorrectionAna art::module SOURCE LightPropagationCorrectionAna_module.cc LIBRARIES ${MODULE_LIBRARIES}) + +build_dictionary(DICTIONARY_LIBRARIES + lardataobj::RecoBase + canvas::canvas +) +add_subdirectory(job) +install_fhicl() +install_source() diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc new file mode 100644 index 000000000..10994c52a --- /dev/null +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc @@ -0,0 +1,150 @@ +//////////////////////////////////////////////////////////////////////// +// Class: SBNDOpT0FinderAna +// Plugin Type: analyzer (art v3_05_01) +// File: SBNDOpT0FinderAna_module.cc +// +// Generated at Mon Oct 12 13:52:37 2020 by Marco Del Tutto using cetskelgen +// from cetlib version v3_10_00. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "sbnobj/SBND/OpFlashTiming/CorrectedOpFlashTiming.h" + +#include "art_root_io/TFileService.h" + +#include "TTree.h" + +#include + +class LightPropagationCorrectionAna; + + +class LightPropagationCorrectionAna : public art::EDAnalyzer { +public: + explicit LightPropagationCorrectionAna(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + LightPropagationCorrectionAna(LightPropagationCorrectionAna const&) = delete; + LightPropagationCorrectionAna(LightPropagationCorrectionAna&&) = delete; + LightPropagationCorrectionAna& operator=(LightPropagationCorrectionAna const&) = delete; + LightPropagationCorrectionAna& operator=(LightPropagationCorrectionAna&&) = delete; + + // Required functions. + void analyze(art::Event const& e) override; + void beginJob() override; +private: + + std::string fLightPropCorrectionLabel; ///< The light propagation correction label producer + + TTree* fTree; ///< The TTree for storing event information + + double fOpFlashT0; + double fUpstreamTime_lightonly; + double fUpstreamTime_tpczcorr; + double fUpstreamTime_propcorr_tpczcorr; + + unsigned int _eventID; + unsigned int _runID; + unsigned int _subrunID; + +}; + + +LightPropagationCorrectionAna::LightPropagationCorrectionAna(fhicl::ParameterSet const& p) + : EDAnalyzer{p} +{ + fLightPropCorrectionLabel = p.get("LightPropCorrectionLabel"); +} + +void LightPropagationCorrectionAna::analyze(art::Event const& e) +{ + fOpFlashT0=-99999.; + fUpstreamTime_lightonly=-99999.; + fUpstreamTime_tpczcorr=-99999.; + fUpstreamTime_propcorr_tpczcorr=-99999.; + + _eventID = -1; + _runID = -1; + _subrunID = -1; + + std::cout << "Run: " << e.id().run() << " Sub: " << e.id().subRun() << " Evt: " << e.id().event() << std::endl; + + _eventID = e.id().event(); + _runID = e.id().run(); + _subrunID = e.id().subRun(); + + // Get all the T0 objects + ::art::Handle> correctedopflash_h; + e.getByLabel(fLightPropCorrectionLabel, correctedopflash_h); + if(!correctedopflash_h.isValid() || correctedopflash_h->empty()) { + mf::LogWarning("fLightPropCorrectionLabel") << "No LightPropCorrectionLabel objects with label " << fLightPropCorrectionLabel << std::endl; + return; + } + + + std::vector> correctedopflash_v; + art::fill_ptr_vector(correctedopflash_v, correctedopflash_h); + + // Get the T0->Slice association + art::FindManyP correctedflash_to_slices(correctedopflash_h, e, fLightPropCorrectionLabel); + art::FindManyP correctedflash_to_flashes(correctedopflash_h, e, fLightPropCorrectionLabel); + + for (size_t i = 0; i < correctedopflash_v.size(); i++) { + + // The T0 object + auto const correctedopflash = correctedopflash_v[i]; + + // The associations from T0 to both slices and flashes + std::vector> slice_v = correctedflash_to_slices.at(correctedopflash.key()); + std::vector> flash_v = correctedflash_to_flashes.at(correctedopflash.key()); + + // One T0 object is always associated to one and only one Slice and Flash + assert(slice_v.size() == 1); + assert(flash_v.size() == 1); + + std::cout << " Corrected flash time is " << correctedopflash->OpFlashT0 << std::endl; + std::cout << " Associated with slice id " << slice_v[0]->ID() << std::endl; + std::cout << "Corrected flash time light only " << correctedopflash->UpstreamTime_lightonly << std::endl; + std::cout << "Corrected flash time tpc z corr " << correctedopflash->UpstreamTime_tpczcorr << std::endl; + std::cout << "Corrected flash time prop corr tpc z corr " << correctedopflash->UpstreamTime_propcorr_tpczcorr << std::endl; + + fOpFlashT0 = correctedopflash->OpFlashT0; + fUpstreamTime_lightonly = correctedopflash->UpstreamTime_lightonly; + fUpstreamTime_tpczcorr = correctedopflash->UpstreamTime_tpczcorr; + fUpstreamTime_propcorr_tpczcorr = correctedopflash->UpstreamTime_propcorr_tpczcorr; + fTree->Fill(); + } +} + +void LightPropagationCorrectionAna::beginJob() +{ + // Implementation of optional member function here. + art::ServiceHandle tfs; + fTree = tfs->make("LightPropagationCorrection", "Light Propagation Correction Tree"); + + fTree->Branch("eventID", &_eventID, "eventID/i"); + fTree->Branch("runID", &_runID, "runID/i"); + fTree->Branch("subrunID", &_subrunID, "subrunID/i"); + + fTree->Branch("fOpFlashT0", &fOpFlashT0, "OpFlashT0/d"); + fTree->Branch("fUpstreamTime_lightonly", &fUpstreamTime_lightonly, "UpstreamTime_lightonly/d"); + fTree->Branch("fUpstreamTime_tpczcorr", &fUpstreamTime_tpczcorr, "UpstreamTime_tpczcorr/d"); + fTree->Branch("fUpstreamTime_propcorr_tpczcorr", &fUpstreamTime_propcorr_tpczcorr, "UpstreamTime_propcorr_tpczcorr/d"); + +} + +DEFINE_ART_MODULE(LightPropagationCorrectionAna) \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc new file mode 100644 index 000000000..6f44c2272 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -0,0 +1,619 @@ +#include "LightPropagationCorrection_module.hh" + +namespace sbnd { + class LightPropagationCorrection; +} + +sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet const& p) + : EDProducer{p}, + fReco2Label( p.get("Reco2Label") ), + fOpT0FinderModuleLabel( p.get("OpT0FinderModuleLabel") ), + fTPCPMTBarycenterFMModuleLabel( p.get("TPCPMTBarycenterFMModuleLabel") ), + fOpFlashLabel_tpc0 ( p.get("OpFlashLabel_tpc0") ), + fOpFlashLabel_tpc1 ( p.get("OpFlashLabel_tpc1") ), + fSpacePointLabel( p.get("SpacePointLabel") ), + fOpHitsModuleLabel( p.get("OpHitsModuleLabel") ), + fSPECTDCLabel( p.get("SPECTDCLabel") ), + fFlashMatchingTool( p.get("FlashMatchingTool") ), + fSaveCorrectionTree( p.get("SaveCorrectionTree") ), + fSaveSPECTDC( p.get("SaveSPECTDC") ), + fSpeedOfLight( p.get("SpeedOfLight") ), + fVGroupVIS( p.get("VGroupVIS") ), + fVGroupVUV( p.get("VGroupVUV") ), + fNuScoreThreshold( p.get("NuScoreThreshold") ), + fFMScoreThreshold( p.get("FMScoreThreshold") ), + fPDFraction( p.get("PDFraction") ), + fPreWindow( p.get("PreWindow") ), + fPostWindow( p.get("PostWindow") ), + fMinHitPE( p.get("MinHitPE") ), + fReadoutDelay( p.get("ReadoutDelay") ), + fDebug( p.get("Debug", false) ) + // + // More initializers here. +{ + // Initialize the TimeCorrectionVector PerChannel + for(size_t i = 0; i < fWireReadout.NOpChannels(); ++i) { + fTimeCorrectionPerChannel.push_back(0.0); // Initialize with zero or any default value + } + + for(unsigned int opch=0; opch()->TPC(); + fDriftDistance = tpc.DriftDistance(); + fKinkDistance = 0.5*fDriftDistance*(1-fVGroupVUV/fVGroupVIS); + fVGroupVUV_I = 1./fVGroupVUV; + fVISLightPropTime = fDriftDistance/fVGroupVIS; + + // Initialize flash algo for both TPCs + auto const flash_algo = p.get("FlashFinderAlgo"); + auto const flash_pset_tpc0 = p.get("AlgoConfig_tpc0"); + auto algo_ptr_tpc0 = ::lightana::FlashAlgoFactory::get().create(flash_algo,flash_algo); + algo_ptr_tpc0->Configure(flash_pset_tpc0); + _mgr_tpc0.SetFlashAlgo(algo_ptr_tpc0); + + auto const flash_pset_tpc1 = p.get("AlgoConfig_tpc1"); + auto algo_ptr_tpc1 = ::lightana::FlashAlgoFactory::get().create(flash_algo,flash_algo); + algo_ptr_tpc1->Configure(flash_pset_tpc1); + _mgr_tpc1.SetFlashAlgo(algo_ptr_tpc1); + + //Initialize flash geo tool + auto const flashgeo_pset = p.get("FlashGeoConfig"); + _flashgeo = art::make_tool(flashgeo_pset); + + //Initialize flash t0 tool + auto const flasht0_pset = p.get("FlashT0Config"); + _flasht0calculator = art::make_tool(flasht0_pset); + + produces< std::vector >(); + produces>(); + produces>(); +} + +void sbnd::LightPropagationCorrection::produce(art::Event & e) +{ + fEvent = e.id().event(); + fRun = e.id().run(); + fSubrun = e.id().subRun(); + + std::unique_ptr< std::vector > correctedOpFlashTimes (new std::vector); + art::PtrMaker make_correctedopflashtime_ptr{e}; + std::unique_ptr< art::Assns> newCorrectedOpFlashTimingSliceAssn (new art::Assns); + std::unique_ptr< art::Assns> newCorrectedOpFlashTimingOpFlashAssn (new art::Assns); + + // --- Read Recob Slice + ::art::Handle> sliceHandle; + e.getByLabel(fReco2Label, sliceHandle); + // Slice to OpT0Finder + //Get the handle for making the assns later on + ::art::Handle> opt0Handle; + e.getByLabel(fOpT0FinderModuleLabel, opt0Handle); + // Slice to TPCPMTBarycenterFM + ::art::Handle> tpcpmtbarycenterfmHandle; + e.getByLabel(fTPCPMTBarycenterFMModuleLabel, tpcpmtbarycenterfmHandle); + + //Read PFPs + ::art::Handle> pfpHandle; + e.getByLabel(fReco2Label, pfpHandle); + //Read OpFlash Handle + art::Handle< std::vector > opflashListHandle_tpc0; + e.getByLabel(fOpFlashLabel_tpc0, opflashListHandle_tpc0); + art::Handle< std::vector > opflashListHandle_tpc1; + e.getByLabel(fOpFlashLabel_tpc1, opflashListHandle_tpc1); + + art::FindManyP slice_opt0finder_assns(sliceHandle, e, fOpT0FinderModuleLabel); + // Slice to TPCPMTBarycenterFM + art::FindManyP slice_tpcpmtbarycentermatching_assns(sliceHandle, e, fTPCPMTBarycenterFMModuleLabel); + // Slice to hits + art::FindManyP slice_hit_assns (sliceHandle, e, fReco2Label); + //Slice to PFParticles association + art::FindManyP slice_pfp_assns (sliceHandle, e, fReco2Label); + //PFP to vertex + art::FindManyP pfp_vertex_assns(pfpHandle, e, fReco2Label); + //PFP to space points + art::FindManyP pfp_sp_assns(pfpHandle, e, fSpacePointLabel); + //OpFlash to OpHit + art::FindManyP flashToOpHitAssns_tpc0(opflashListHandle_tpc0, e, fOpFlashLabel_tpc0); + art::FindManyP flashToOpHitAssns_tpc1(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1); + // PFP Metadata + art::FindManyP pfp_to_metadata(pfpHandle, e, fReco2Label); + + + // --- Store candidate slices + std::vector< art::Ptr > sliceVect; + art::fill_ptr_vector(sliceVect, sliceHandle); + + //Vector for recob PFParticles + std::vector> pfpVect; + + // --- Get the candidate slices + for(size_t ix=0; ix &pfp : pfpVect){ + if(pfp->IsPrimary() &&( std::abs(pfp->PdgCode())==12 || std::abs(pfp->PdgCode())==14 ) ){ + const std::vector> pfpMetaVec = pfp_to_metadata.at(pfp->Self()); + for (auto const pfpMeta : pfpMetaVec) { + larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap(); + _fNuScore = propertiesMap.at("NuScore"); + if(_fNuScore>_sliceMaxNuScore) _sliceMaxNuScore = _fNuScore; + } + } + + std::vector< art::Ptr > vertexVec = pfp_vertex_assns.at(pfp.key()); + for(const art::Ptr &ver : vertexVec){ + geo::Point_t xyz_vertex = ver->position(); + fRecoVx= xyz_vertex.X(); + fRecoVy= xyz_vertex.Y(); + fRecoVz= xyz_vertex.Z(); + } + + // If correct light propagation time + + + // Get the SP associated to the PFP and then get the hits associated to the SP. ---> Hits associated to the PFP + //Get the spacepoints associated to the PFParticle + std::vector> PFPSpacePointsVect = pfp_sp_assns.at(pfp.key()); + //Get the SP Hit assns + art::Handle> eventSpacePoints; + std::vector> eventSpacePointsVect; + e.getByLabel(fSpacePointLabel, eventSpacePoints); + art::fill_ptr_vector(eventSpacePointsVect, eventSpacePoints); + art::FindManyP SPToHitAssoc (eventSpacePointsVect, e, fSpacePointLabel); + for (const art::Ptr &SP: PFPSpacePointsVect){ + std::vector> SPHit = SPToHitAssoc.at(SP.key()); + if (SPHit.at(0)->WireID().Plane==2){ + fSpacePointX.push_back(SP->position().X()); + fSpacePointY.push_back(SP->position().Y()); + fSpacePointZ.push_back(SP->position().Z()); + fSpacePointIntegral.push_back(SPHit.at(0)->Integral()); + + //Fill Bayrcenter Position + if(SP->position().X() < 0){ + fChargeWeightX[0] += SP->position().X() * SPHit.at(0)->Integral(); + fChargeWeightY[0] += SP->position().Y() * SPHit.at(0)->Integral(); + fChargeWeightZ[0] += SP->position().Z() * SPHit.at(0)->Integral(); + fChargeTotalWeight[0] += SPHit.at(0)->Integral(); + } + else{ + fChargeWeightX[1] += SP->position().X() * SPHit.at(0)->Integral(); + fChargeWeightY[1] += SP->position().Y() * SPHit.at(0)->Integral(); + fChargeWeightZ[1] += SP->position().Z() * SPHit.at(0)->Integral(); + fChargeTotalWeight[1] += SPHit.at(0)->Integral(); + } + } + } + } + //Fill TPC 0 information + fChargeBarycenterX[0] = fChargeWeightX[0]/fChargeTotalWeight[0]; + fChargeBarycenterY[0] = fChargeWeightY[0]/fChargeTotalWeight[0]; + fChargeBarycenterZ[0] = fChargeWeightZ[0]/fChargeTotalWeight[0]; + //Fill TPC 1 information + fChargeBarycenterX[1] = fChargeWeightX[1]/fChargeTotalWeight[1]; + fChargeBarycenterY[1] = fChargeWeightY[1]/fChargeTotalWeight[1]; + fChargeBarycenterZ[1] = fChargeWeightZ[1]/fChargeTotalWeight[1]; + + if(_sliceMaxNuScoreGetPropagationTimeCorrectionPerChannel(); + + // Get the SPECTDC product required to go to the RWM reference frame + if(fSaveSPECTDC) + { + art::Handle> tdcHandle; + e.getByLabel(fSPECTDCLabel, tdcHandle); + if (!tdcHandle.isValid() || tdcHandle->size() == 0){ + std::cout << "No SPECTDC products found. Skip this event." << std::endl; + return; + } + else{ + const std::vector tdc_v(*tdcHandle); + for (size_t i=0; i> flashFM; + if(fFlashMatchingTool == "OpT0Finder" ){ + const std::vector< art::Ptr > slcOpT0Finder = slice_opt0finder_assns.at( slice.key() ); + if(slcOpT0Finder.size() == 0) continue; + size_t OpT0Idx = HighestOpT0ScoreIdx(slcOpT0Finder); + // Get the flash OpT0 association + art::FindManyP opflash_opt0finder_assns(opt0Handle, e, fOpT0FinderModuleLabel); + flashFM = opflash_opt0finder_assns.at( slcOpT0Finder[OpT0Idx].key() ); + if(flashFM.size() > 1){ + throw art::Exception(art::errors::LogicError) << "There are multiple OpFlash objects associated to the same OpT0Finder object. This is not expected."; + } + _fFMScore = slcOpT0Finder[OpT0Idx]->score; + if(_fFMScore < fFMScoreThreshold){ + ResetSliceInfo(); + continue; + } + } + else if(fFlashMatchingTool == "BarycenterFM") + { + const std::vector< art::Ptr > slcTPCPMTBarycenter = slice_tpcpmtbarycentermatching_assns.at( slice.key() ); + if(slcTPCPMTBarycenter.size() == 0) continue; + art::FindManyP opflash_tpcpmtbarycenter_assns(tpcpmtbarycenterfmHandle, e, fTPCPMTBarycenterFMModuleLabel); + size_t BFMIdx = HighestBFMScoreIdx(slcTPCPMTBarycenter); + flashFM = opflash_tpcpmtbarycenter_assns.at( slcTPCPMTBarycenter[BFMIdx].key() ); + if(flashFM.size() > 1){ + throw art::Exception(art::errors::LogicError) << "There are multiple OpFlash objects associated to the same TPCPMTBarycenterFM object. This is not expected."; + } + _fFMScore = slcTPCPMTBarycenter[BFMIdx]->score; + if(_fFMScore < fFMScoreThreshold){ + ResetSliceInfo(); + continue; + } + } + else throw art::Exception(art::errors::LogicError) << " Flash matching tool " << fFlashMatchingTool << " not supported ." << std::endl; + + + // Get the ophits associated to the flash + std::vector> ophitlist; + if(flashFM[0]->XCenter()<0) + { + ophitlist = flashToOpHitAssns_tpc0.at(flashFM[0].key()); + _mgr = _mgr_tpc0; // Use the TPC 0 flash finder manager + } + else + { + ophitlist = flashToOpHitAssns_tpc1.at(flashFM[0].key()); + _mgr = _mgr_tpc1; // Use the TPC 1 flash finder manager + } + + std::vector newOpHitList; + std::vector oldOpHitList; + for(const auto& ophit : ophitlist) { + oldOpHitList.push_back(*ophit); + } + // Get the list of the corrected OpHits + this->CorrectOpHitTime(ophitlist, newOpHitList); + // Create the list of ophit lite to be used in the flash finder + ::lightana::LiteOpHitArray_t ophits; + this->FillLiteOpHit(newOpHitList, ophits); + // Create the flash manager + auto const flash_v = _mgr.RecoFlash(ophits); + double originalFlashTime = flashFM[0]->Time(); + double newFlashTime = 0.0; + for(const auto& lflash : flash_v) { + // Get Flash Barycenter + double Ycenter, Zcenter, Ywidth, Zwidth; + _flashgeo->GetFlashLocation(lflash.channel_pe, Ycenter, Zcenter, Ywidth, Zwidth); + // Get flasht0 + double flasht0 = lflash.time; + // Refine t0 calculation + flasht0 = _flasht0calculator->GetFlashT0(lflash.time, GetAssociatedLiteHits(lflash, ophits)); + // Subtract readout ReadoutDelay + flasht0 = flasht0 - fReadoutDelay; + recob::OpFlash flash(flasht0, lflash.time_err, flasht0, + ( flasht0) / 1600., lflash.channel_pe, + 0, 0, 1, // this are just default values + 100., -1., Ycenter, Ywidth, Zcenter, Zwidth); + newFlashTime = flasht0; + sbnd::OpFlashTiming::CorrectedOpFlashTiming correctedOpFlashTiming; + correctedOpFlashTiming.OpFlashT0 = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000; + correctedOpFlashTiming.UpstreamTime_lightonly = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (Zcenter/fSpeedOfLight)/1000; + correctedOpFlashTiming.UpstreamTime_tpczcorr = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (fRecoVz/fSpeedOfLight)/1000; + correctedOpFlashTiming.UpstreamTime_propcorr_tpczcorr = newFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (fRecoVz/fSpeedOfLight)/1000; + correctedOpFlashTiming.OpT0Score = _fFMScore; + correctedOpFlashTiming.SliceNuScore = _sliceMaxNuScore; + correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming)); + } + + art::Ptr newCorrectedOpFlashTimingPtr = make_correctedopflashtime_ptr(correctedOpFlashTimes->size()-1); + newCorrectedOpFlashTimingSliceAssn->addSingle(slice, newCorrectedOpFlashTimingPtr); + newCorrectedOpFlashTimingOpFlashAssn->addSingle(flashFM[0], newCorrectedOpFlashTimingPtr); + if(fSaveCorrectionTree){ + this->FillCorrectionTree(newFlashTime, *flashFM[0], oldOpHitList, newOpHitList); + } + } + + fTree->Fill(); + ResetEventVars(); + e.put(std::move(correctedOpFlashTimes)); + e.put(std::move(newCorrectedOpFlashTimingSliceAssn)); + e.put(std::move(newCorrectedOpFlashTimingOpFlashAssn)); + return; +} + +void sbnd::LightPropagationCorrection::beginJob() +{ + if(fSaveCorrectionTree) + { + fTree = tfs->make("PMTWaveformFilteranalyzer", "PMT Waveform Filter Analyzer Tree"); + fTree->Branch("eventID", &fEvent, "eventID/i"); + fTree->Branch("runID", &fRun, "runID/i"); + fTree->Branch("subrunID", &fSubrun, "subrunID/i"); + fTree->Branch("RWMTime", &fRWMTime); + fTree->Branch("EventTriggerTime", &fEventTriggerTime); + fTree->Branch("NuScore", &fNuScore); + fTree->Branch("FMScore", &fFMScore); + fTree->Branch("OpFlashTimeOld", &fOpFlashTimeOld); + fTree->Branch("OpFlashTimeNew", &fOpFlashTimeNew); + fTree->Branch("OpFlashXCenter", &fOpFlashXCenter); + fTree->Branch("OpFlashYCenter", &fOpFlashYCenter); + fTree->Branch("OpFlashZCenter", &fOpFlashZCenter); + fTree->Branch("OpFlashPE", &fOpFlashPE); + fTree->Branch("SliceVx", &fSliceVx); + fTree->Branch("SliceVy", &fSliceVy); + fTree->Branch("SliceVz", &fSliceVz); + fTree->Branch("SliceSPX", &fSliceSPX); + fTree->Branch("SliceSPY", &fSliceSPY); + fTree->Branch("SliceSPZ", &fSliceSPZ); + fTree->Branch("OpHitOldTime", &fOpHitOldTime); + fTree->Branch("OpHitNewTime", &fOpHitNewTime); + fTree->Branch("OpHitPE", &fOpHitPE); + fTree->Branch("OpHitOpCh", &fOpHitOpCh); + } +} + +void sbnd::LightPropagationCorrection::endJob() +{ + +} + +void sbnd::LightPropagationCorrection::ResetEventVars() +{ + if(fSaveCorrectionTree) + { + fEvent = 0; + fRun = 0; + fSubrun = 0; + _fNuScore = 0.0; + fRWMTime=-99999.; + fEventTriggerTime=-99999.; + fNuScore.clear(); + fFMScore.clear(); + fOpFlashTimeOld.clear(); + fOpFlashTimeNew.clear(); + fOpFlashXCenter.clear(); + fOpFlashYCenter.clear(); + fOpFlashZCenter.clear(); + fOpHitOldTime.clear(); + fOpHitNewTime.clear(); + fOpHitPE.clear(); + fOpHitOpCh.clear(); + fOpFlashPE.clear(); + fSliceVx.clear(); + fSliceVy.clear(); + fSliceVz.clear(); + fSliceSPX.clear(); + fSliceSPY.clear(); + fSliceSPZ.clear(); + } +} + +size_t sbnd::LightPropagationCorrection::HighestOpT0ScoreIdx(const std::vector< art::Ptr > slcFM) +{ + // Gets the idx of the OpT0 object with the highest score + double highestOpT0Score = -99999.0; // Initialize to a negative value + size_t highestIdx = 0; + for(size_t jx=0; jxscore > highestOpT0Score){ + highestOpT0Score = slcFM[jx]->score; + highestIdx = jx; + } + } + return highestIdx; +} + +size_t sbnd::LightPropagationCorrection::HighestBFMScoreIdx(const std::vector< art::Ptr > slcFM) +{ + // Gets the idx of the OpT0 object with the highest score + double highesBFM0Score = -99999.0; // Initialize to a negative value + size_t highestIdx = 0; + for(size_t jx=0; jxscore > highesBFM0Score){ + highesBFM0Score = slcFM[jx]->score; + highestIdx = jx; + } + } + return highestIdx; +} + +void sbnd::LightPropagationCorrection::ResetSliceInfo() +{ + _fNuScore=0.0; + fRecoVx = 0.0; + fRecoVy = 0.0; + fRecoVz = 0.0; + fSpacePointX.clear(); + fSpacePointY.clear(); + fSpacePointZ.clear(); + fSpacePointIntegral.clear(); + fTimeCorrectionPerChannel.resize(312, 0.0); // Reset the time correction vector for each channel + fChargeBarycenterX.assign(2, 0.0); + fChargeBarycenterY.assign(2, 0.0); + fChargeBarycenterZ.assign(2, 0.0); + fChargeWeightX.assign(2, 0.0); + fChargeWeightY.assign(2, 0.0); + fChargeWeightZ.assign(2, 0.0); + fChargeTotalWeight.assign(2, 0.0); +} + +void sbnd::LightPropagationCorrection::GetPropagationTimeCorrectionPerChannel() +{ + // Implementation + for(size_t opdet = 0; opdet < fOpDetID.size(); ++opdet) { + double _opDetX = fOpDetX[opdet]; + double _opDetY = fOpDetY[opdet]; + double _opDetZ = fOpDetZ[opdet]; + float minPropTime = 999999999.; + + for(size_t sp=0; sp fKinkDistance) + return (fDriftDistance-std::abs(drift)) * fVGroupVUV_I ; + else + return std::abs(drift) * fVGroupVUV_I + fVISLightPropTime; +} + + +double sbnd::LightPropagationCorrection::GetFlashT0(double flash_time, std::vector ophitlist){ + + std::vector< std::pair > selected_hits; + double pe_sum = 0; + + // fill vector with selected hits in the specified window + for(auto const& hit : ophitlist) { + if( hit.PeakTime()flash_time-fPreWindow && hit.PE()>fMinHitPE){ + selected_hits.push_back( std::make_pair(hit.PE(), hit.PeakTime())); + pe_sum += hit.PE(); + } + } + + if(pe_sum>0){ + // sort vector by number of #PE (ascending order) + std::sort( selected_hits.begin(), selected_hits.end(), std::greater< std::pair >() ); + + double flasht0_mean=0., pe_count=0.; + int nophits=0; + + // loop over selected ophits + for (size_t ix=0; ixfPDFraction ) break; + } + return flasht0_mean/nophits; + } + else + return flash_time; +} + + +void sbnd::LightPropagationCorrection::CorrectOpHitTime(std::vector> OldOpHitList, std::vector & newOpHitList) +{ + int _nophits = 0; + _nophits += OldOpHitList.size(); + for (int i = 0; i < _nophits; ++i) { + int opCh = OldOpHitList.at(i)->OpChannel(); + double channelCorrection = fTimeCorrectionPerChannel[opCh]; + double newPeakTime = OldOpHitList.at(i)->StartTime()+OldOpHitList.at(i)->RiseTime() + channelCorrection/1000; + double newPeakTimeAbs = OldOpHitList.at(i)->PeakTimeAbs()+ channelCorrection/1000; + double newStartTime = OldOpHitList.at(i)->StartTime() + channelCorrection/1000; + double riseTime = OldOpHitList.at(i)->RiseTime(); + unsigned int frame = OldOpHitList.at(i)->Frame(); + double width = OldOpHitList.at(i)->Width(); + double area = OldOpHitList.at(i)->Area(); + double amplitude = OldOpHitList.at(i)->Amplitude(); + double pe = OldOpHitList.at(i)->PE(); + recob::OpHit newOpHit = recob::OpHit(opCh, newPeakTime, newPeakTimeAbs, newStartTime, riseTime, frame, width, area, amplitude, pe, 0.0); + newOpHitList.push_back(newOpHit); + } +} + + +void sbnd::LightPropagationCorrection::FillLiteOpHit(std::vector const& OpHitList, std::vector<::lightana::LiteOpHit_t>& LiteOpHitList) +{ + for(auto const& oph : OpHitList) { + ::lightana::LiteOpHit_t loph; + loph.peak_time = oph.StartTime()+oph.RiseTime(); + loph.pe = oph.PE(); + loph.channel = oph.OpChannel(); + LiteOpHitList.emplace_back(std::move(loph)); + } +} + + +void sbnd::LightPropagationCorrection::FillCorrectionTree(double & newFlashTime, recob::OpFlash const& flash, std::vector const& oldOpHitList, std::vector const& newOpHitList){ + + fSliceSPX.push_back({}); + fSliceSPY.push_back({}); + fSliceSPZ.push_back({}); + fOpHitOldTime.push_back({}); + fOpHitNewTime.push_back({}); + fOpHitPE.push_back({}); + fOpHitOpCh.push_back({}); + + if(fDebug) + { + for(size_t i=0; i +#include + +// Services +#include "larsim/MCCheater/BackTrackerService.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "larsim/Utils/TruthMatchUtils.h" +#include "lardataalg/DetectorInfo/DetectorClocksData.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "larcore/Geometry/Geometry.h" +#include "larcore/Geometry/WireReadout.h" +#include "lardata/Utilities/AssociationUtil.h" + + +// G4 includes +#include "nusimdata/SimulationBase/MCParticle.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "lardataobj/Simulation/SimPhotons.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" + +// Reco includes +// PDS +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "lardataobj/RawData/OpDetWaveform.h" +// TPC +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/SpacePoint.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCNeutrino.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "lardataobj/RawData/RawDigit.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/Shower.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/SpacePoint.h" +#include "lardataobj/RecoBase/Cluster.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RawData/RawDigit.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/RecoBase/Vertex.h" +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/AnalysisBase/Calorimetry.h" +#include "lardataobj/AnalysisBase/ParticleID.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "lardataobj/RecoBase/PFParticleMetadata.h" + +// CRT +#include "sbnobj/SBND/CRT/CRTTrack.hh" +#include "sbnobj/SBND/OpFlashTiming/CorrectedOpFlashTiming.h" + +// Cosmic rejection includes +#include "sbnobj/Common/Reco/OpT0FinderResult.h" +#include "sbnobj/Common/Reco/SimpleFlashMatchVars.h" +#include "sbnobj/Common/Reco/CRUMBSResult.h" +#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" +#include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "lardataobj/AnalysisBase/T0.h" + +// Geometry and mapping +#include "larcore/Geometry/WireReadout.h" +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" +#include +#include + +// Flash finder utilities +#include "sbndcode/OpDetReco/OpFlash/FlashFinder/FlashFinderManager.h" +#include "sbndcode/OpDetReco/OpFlash/FlashFinder/FlashFinderFMWKInterface.h" +#include "sbndcode/OpDetReco/OpFlash/FlashFinder/PECalib.h" +#include "sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoBase.hh" +#include "sbndcode/OpDetReco/OpFlash/FlashTools/FlashT0Base.hh" +#include "sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorBase.hh" + +#define fXFidCut1 1.5 +#define fXFidCut2 190 +#define fYFidCut 190 +#define fZFidCut1 10 +#define fZFidCut2 490 + +#define xdet_size 1000 +#define ydet_size 1000 +#define zmindet_size -500 +#define zmaxdet_size 1800 + +#define fDefaulNeutrinoID 99999 + + +namespace sbnd { + class LightPropagationCorrection; + +} + +class sbnd::LightPropagationCorrection : public art::EDProducer { +public: + explicit LightPropagationCorrection(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + LightPropagationCorrection(LightPropagationCorrection const&) = delete; + LightPropagationCorrection(LightPropagationCorrection&&) = delete; + LightPropagationCorrection& operator=(LightPropagationCorrection const&) = delete; + LightPropagationCorrection& operator=(LightPropagationCorrection&&) = delete; + + +private: + + // Required functions. + void produce(art::Event & e) override; + void endJob() override; + void beginJob() override; + + // Selected optional functions. + void ResetEventVars(); + void ResetSliceInfo(); + size_t HighestOpT0ScoreIdx(const std::vector< art::Ptr >); + size_t HighestBFMScoreIdx(const std::vector< art::Ptr > ); + + void GetPropagationTimeCorrectionPerChannel(); + void CorrectOpHitTime(std::vector> , std::vector & ); + void FillLiteOpHit(std::vector const& , std::vector<::lightana::LiteOpHit_t>& ); + double GetPropagationTime(double ); + double GetFlashT0(double , std::vector ); + void FillCorrectionTree(double & , recob::OpFlash const& , std::vector const& , std::vector const& ); + ::lightana::LiteOpHitArray_t GetAssociatedLiteHits(::lightana::LiteOpFlash_t , ::lightana::LiteOpHitArray_t ); + + + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); + + //Flash finder manager + ::lightana::FlashFinderManager _mgr; + ::lightana::FlashFinderManager _mgr_tpc0; + ::lightana::FlashFinderManager _mgr_tpc1; + + // Tool for calculating the OpFlash Y and Z centers + std::unique_ptr _flashgeo; + + // Tool for calculating the OpFlash t0 + std::unique_ptr _flasht0calculator; + + //Vector for PMT position + std::vector fOpDetID; + std::vector fOpDetX; + std::vector fOpDetY; + std::vector fOpDetZ; + + std::string fReco2Label; + std::string fOpT0FinderModuleLabel; + std::string fTPCPMTBarycenterFMModuleLabel; + std::string fOpFlashLabel_tpc0; + std::string fOpFlashLabel_tpc1; + std::string fSpacePointLabel; + std::string fOpHitsModuleLabel; + std::string fOpFlashNewLabel; + std::string fSPECTDCLabel; + std::string fFlashMatchingTool; + + bool fSaveCorrectionTree; + bool fSaveDebugTree; + bool fSaveSPECTDC; + + std::vector fTimeCorrectionPerChannel; + double fRecoVx = 0.0; + double fRecoVy = 0.0; + double fRecoVz = 0.0; + + //Space Point Variables + std::vector fSpacePointX; + std::vector fSpacePointY; + std::vector fSpacePointZ; + std::vector fSpacePointIntegral; + + //Charge Barycenter + std::vector fChargeBarycenterX{0.,0.}; + std::vector fChargeBarycenterY{0.,0.}; + std::vector fChargeBarycenterZ{0.,0.}; + std::vector fChargeWeightX{0.,0.}; + std::vector fChargeWeightY{0.,0.}; + std::vector fChargeWeightZ{0.,0.}; + std::vector fChargeTotalWeight{0.,0.}; + + double fDriftDistance; // Total Drift Distance + double fSpeedOfLight; // Speed of light in mm/ns + double fVISLightPropTime; + double fKinkDistance; + double fVGroupVUV_I; + double fVGroupVIS; + double fVGroupVUV; // Speed of light in vacuum in mm/ns + + double fNuScoreThreshold; + double fFMScoreThreshold; + + double fPDFraction; + double fPreWindow; + double fPostWindow; + double fMinHitPE; + double fReadoutDelay; + + bool fDebug; + + art::ServiceHandle tfs; + TTree *fTree; + + int fEvent; + int fRun; + int fSubrun; + double _fNuScore; + double _fFMScore; + double fEventTriggerTime=-999999.; + double fRWMTime=-999999.; + std::vector fNuScore; + std::vector fFMScore; + std::vector fOpFlashTimeOld; + std::vector fOpFlashTimeNew; + std::vector fOpFlashXCenter; + std::vector fOpFlashYCenter; + std::vector fOpFlashZCenter; + std::vector fOpFlashPE; + std::vector fSliceVx; + std::vector fSliceVy; + std::vector fSliceVz; + std::vector> fSliceSPX; + std::vector> fSliceSPY; + std::vector> fSliceSPZ; + std::vector> fOpHitOldTime; + std::vector> fOpHitNewTime; + std::vector> fOpHitPE; + std::vector> fOpHitOpCh; +}; + + +DEFINE_ART_MODULE(sbnd::LightPropagationCorrection) diff --git a/sbndcode/LightPropagationCorrection/classes.h b/sbndcode/LightPropagationCorrection/classes.h new file mode 100644 index 000000000..acc2700b8 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/classes.h @@ -0,0 +1,15 @@ +//File: classes.h +//Brief: Include directives needed to generate header(s) for the raw::pmt::sbndpmt data product. +//Author: Lynn Tung + +//ART includes +#include "canvas/Persistency/Common/Wrapper.h" +#include "canvas/Persistency/Common/Assns.h" + +//local includes +#include "sbndcode/Timing/SBNDRawTimingObj.h" + +//larsoft includes +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardataobj/RecoBase/Slice.h" diff --git a/sbndcode/LightPropagationCorrection/classes_def.xml b/sbndcode/LightPropagationCorrection/classes_def.xml new file mode 100644 index 000000000..3954e82d3 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/classes_def.xml @@ -0,0 +1,12 @@ + + + + + + + + diff --git a/sbndcode/LightPropagationCorrection/job/CMakeLists.txt b/sbndcode/LightPropagationCorrection/job/CMakeLists.txt new file mode 100644 index 000000000..eca2a52d6 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/job/CMakeLists.txt @@ -0,0 +1 @@ +install_fhicl() \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl new file mode 100644 index 000000000..218248c50 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl @@ -0,0 +1,44 @@ +#include "services_sbnd.fcl" +#include "sbnd_lightpropagationcorrection_config.fcl" + +process_name: LIGHTPROPAGATIONCORR # The process name must NOT contain any underscores + +source: +{ + module_type: RootInput # Telling art we want a ROOT input + maxEvents: -1 +} + +services: +{ + TFileService: { fileName: "lightpropagationcorrection.root" } + @table::sbnd_services +} + +outputs: +{ + out1: + { + fileName: "%ifb_lightpropcorr.root" + module_type: RootOutput + dataTier: "reco" + #outputCommands: [ + # "keep *_*_*_*" + # , "drop *_pmtdecoder_*_*" + # ] + } +} + +physics: +{ + producers: + { + lightpropcorrection: @local::LightPropagationCorrection + } + + prod: [ lightpropcorrection ] + stream1: [ out1 ] + + trigger_paths: [ prod ] + end_paths: [ stream1 ] +} diff --git a/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl new file mode 100644 index 000000000..52992e020 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl @@ -0,0 +1,33 @@ +#include "services_sbnd.fcl" +#include "simulationservices_sbnd.fcl" +#include "messages_sbnd.fcl" +#include "sam_sbnd.fcl" +#include "sbnd_lightpropagationcorrectionana_config.fcl" + +process_name: CorrectedOpFlashAna + +services: +{ + TFileService: { fileName: "correctedopflashana_tree.root" } + @table::sbnd_basic_services + @table::sbnd_simulation_services + ParticleInventoryService: @local::sbnd_particleinventoryservice + BackTrackerService: @local::sbnd_backtrackerservice + ParticleInventoryService: @local::standard_particleinventoryservice +} + +source: +{ + module_type: RootInput + maxEvents: -1 # Number of events to create +} + +physics: +{ + analyzers: + { + lightpropcorrectionanatree: @local::lightpropagationcorrection_ana + } + ana: [lightpropcorrectionanatree] + end_paths: [ana] +} \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl new file mode 100644 index 000000000..04788fb75 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl @@ -0,0 +1,49 @@ +#include "opticalsimparameterisations_sbnd.fcl" +#include "sbnd_flashalgo.fcl" +#include "sbnd_flasht0algo.fcl" +#include "sbnd_flashgeoalgo.fcl" +#include "sbnd_flashfinder.fcl" + +BEGIN_PROLOG + +LightPropagationCorrection: +{ + module_type: "LightPropagationCorrection" + OpFlashLabel_tpc0: "opflashtpc0" + OpFlashLabel_tpc1: "opflashtpc1" + + Reco2Label: "pandora" + OpHitsModuleLabel: "ophitpmt" + SpacePointLabel: "pandora" + OpT0FinderModuleLabel: "opt0finder" + TPCPMTBarycenterFMModuleLabel: "tpcpmtbarycentermatching" + SPECTDCLabel: "tdcdecoder" + + VGroupVUV: @local::sbnd_vuv_timing_parameterization.vuv_vgroup_mean + VGroupVIS: @local::sbnd_vis_timing_parameterization.vis_vmean + SpeedOfLight: 29.9 + + FlashT0Config: @local::FlashT0SelectedChannels + MinHitPE: @local::FlashT0SelectedChannels.MinHitPE + PreWindow: @local::FlashT0SelectedChannels.PreWindow + PostWindow: @local::FlashT0SelectedChannels.PostWindow + PDFraction: @local::FlashT0SelectedChannels.PDFraction + FlashGeoConfig: @local::FlashGeoThreshold + AlgoConfig_tpc0: @local::SimpleFlashStandard + AlgoConfig_tpc1: @local::SimpleFlashStandard + FlashFinderAlgo: @local::SBNDSimpleFlash.FlashFinderAlgo + ReadoutDelay: 1.35e-1 + + SaveSPECTDC: true + SaveCorrectionTree: true + + NuScoreThreshold: 0.0 + FMScoreThreshold: 0.0 + + FlashMatchingTool: "BarycenterFM" +} + +LightPropagationCorrection.AlgoConfig_tpc0.TPC: 0 +LightPropagationCorrection.AlgoConfig_tpc1.TPC: 1 + +END_PROLOG \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl new file mode 100644 index 000000000..2d0ca7013 --- /dev/null +++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl @@ -0,0 +1,10 @@ + +BEGIN_PROLOG + +lightpropagationcorrection_ana: +{ + module_type: "LightPropagationCorrectionAna" + LightPropCorrectionLabel: "lightpropcorrection" +} + +END_PROLOG \ No newline at end of file diff --git a/sbndcode/OpDetAnalyzer/PDSAnalyzer/CMakeLists.txt b/sbndcode/OpDetAnalyzer/PDSAnalyzer/CMakeLists.txt index 5ae5a3220..f8e3ac3a8 100644 --- a/sbndcode/OpDetAnalyzer/PDSAnalyzer/CMakeLists.txt +++ b/sbndcode/OpDetAnalyzer/PDSAnalyzer/CMakeLists.txt @@ -23,7 +23,14 @@ set(MODULE_LIBRARIES nusimdata::SimulationBase nug4::ParticleNavigation + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_SBND + sbndaq_artdaq_core::sbndaq-artdaq-core_Obj_SBND + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + sbnobj::Common_Reco + sbnobj::SBND_CRT + sbnobj::SBND_Timing sbndcode_OpDetSim ) diff --git a/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.cc b/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.cc index 95abd309e..0461378b5 100644 --- a/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.cc +++ b/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.cc @@ -12,8 +12,14 @@ opdet::SBNDPDSAnalyzer::SBNDPDSAnalyzer(fhicl::ParameterSet const& p) fSaveRawWaveforms( p.get("SaveRawWaveforms") ), fSaveDeconvolvedWaveforms( p.get("SaveDeconvolvedWaveforms") ), fSaveOpHits( p.get("SaveOpHits") ), + fSaveOnlyFlashHits( p.get("SaveOnlyFlashHits") ), fSaveOpFlashes( p.get("SaveOpFlashes") ), + fSaveCRT( p.get("SaveCRT") ), + fSaveSPECTDC( p.get("SaveSPECTDC") ), + fSaveOnlyCRTPDSMatch( p.get("SaveOnlyCRTPDSMatch") ), + fSaveOnlyAVTracks( p.get("SaveOnlyAVTracks") ), fSaveCosmicId( p.get("SaveCosmicId") ), + fSavePEFlavourPerFlash( p.get("SavePEFlavourPerFlash") ), fVerbosity( p.get("Verbosity") ), fMakePerTrackTree( p.get("MakePerTrackTree") ), fMakePDSGeoTree( p.get("MakePDSGeoTree") ), @@ -30,6 +36,8 @@ opdet::SBNDPDSAnalyzer::SBNDPDSAnalyzer(fhicl::ParameterSet const& p) fDeconvolvedWaveformsModuleLabel( p.get< std::string >("DeconvolvedWaveformsModuleLabel") ), fOpHitsModuleLabel( p.get>("OpHitsModuleLabel") ), fOpFlashesModuleLabel( p.get>("OpFlashesModuleLabel") ), + fCRTTrackModuleLabel( p.get("CRTTrackModuleLabel") ), + fSPECTDCLabel( p.get< std::string >("SPECTDCLabel") ), fHitsLabel( p.get("HitsLabel") ), fReco2Label( p.get("Reco2Label") ), fCosmicIdModuleLabel( p.get< std::string >("CosmicIdModuleLabel") ), @@ -38,7 +46,9 @@ opdet::SBNDPDSAnalyzer::SBNDPDSAnalyzer(fhicl::ParameterSet const& p) fG4BufferBoxX( p.get>("G4BufferBoxX") ), fG4BufferBoxY( p.get>("G4BufferBoxY") ), fG4BufferBoxZ( p.get>("G4BufferBoxZ") ), - fG4BeamWindow( p.get>("G4BeamWindow") ) + fG4BeamWindow( p.get>("G4BeamWindow") ), + fOpFlashBeamWindow( p.get>("OpFlashBeamWindow") ), + fCRTSaveWindow( p.get>("CRTSaveWindow") ) { // MakePerTrackTree mode requires UseSimPhotonsLite false and SaveMCParticles true @@ -166,6 +176,10 @@ void opdet::SBNDPDSAnalyzer::beginJob() fTree->Branch("flash_zerr", "std::vector", &_flash_zerr); fTree->Branch("flash_x","std::vector", &_flash_x); fTree->Branch("flash_xerr", "std::vector", &_flash_xerr); + fTree->Branch("flash_pe_co", "std::vector", &_flash_pe_co); + fTree->Branch("flash_pe_co", "std::vector", &_flash_pe_co); + fTree->Branch("flash_pe_unco", "std::vector", &_flash_pe_unco); + fTree->Branch("flash_pmt_ratio", "std::vector", &_flash_pmt_ratio); fTree->Branch("flash_ophit_time", "std::vector>", &_flash_ophit_time); fTree->Branch("flash_ophit_risetime", "std::vector>", &_flash_ophit_risetime); fTree->Branch("flash_ophit_starttime", "std::vector>", &_flash_ophit_starttime); @@ -176,6 +190,40 @@ void opdet::SBNDPDSAnalyzer::beginJob() fTree->Branch("flash_ophit_ch", "std::vector>", &_flash_ophit_ch); } + if(fSaveCRT) + { + fTree->Branch("tr_entry_x", "std::vector", &_tr_entry_x); + fTree->Branch("tr_entry_y", "std::vector", &_tr_entry_y); + fTree->Branch("tr_entry_z", "std::vector", &_tr_entry_z); + fTree->Branch("tr_exit_x", "std::vector", &_tr_exit_x); + fTree->Branch("tr_exit_y", "std::vector", &_tr_exit_y); + fTree->Branch("tr_exit_z", "std::vector", &_tr_exit_z); + fTree->Branch("tr_start_x", "std::vector", &_tr_start_x); + fTree->Branch("tr_start_y", "std::vector", &_tr_start_y); + fTree->Branch("tr_start_z", "std::vector", &_tr_start_z); + fTree->Branch("tr_end_x", "std::vector", &_tr_end_x); + fTree->Branch("tr_end_y", "std::vector", &_tr_end_y); + fTree->Branch("tr_end_z", "std::vector", &_tr_end_z); + fTree->Branch("tr_dir_x", "std::vector", &_tr_dir_x); + fTree->Branch("tr_dir_y", "std::vector", &_tr_dir_y); + fTree->Branch("tr_dir_z", "std::vector", &_tr_dir_z); + fTree->Branch("tr_ts0", "std::vector", &_tr_ts0); + fTree->Branch("tr_ets0", "std::vector", &_tr_ets0); + fTree->Branch("tr_ets1", "std::vector", &_tr_ets1); + fTree->Branch("tr_ets1", "std::vector", &_tr_ets1); + fTree->Branch("tr_pe", "std::vector", &_tr_pe); + fTree->Branch("tr_length", "std::vector", &_tr_length); + fTree->Branch("tr_tof", "std::vector", &_tr_tof); + fTree->Branch("tr_theta", "std::vector", &_tr_theta); + fTree->Branch("tr_phi", "std::vector", &_tr_phi); + fTree->Branch("tr_triple", "std::vector", &_tr_triple); + } + + if(fSaveSPECTDC){ + fTree->Branch("event_trigger_time", &event_trigger_time); + fTree->Branch("rwm_time", &rwm_time); + } + // Cosmic ID if(fSaveCosmicId){ fTree->Branch("CRUMBSScore","std::vector", &_CRUMBSScore); @@ -242,6 +290,10 @@ void opdet::SBNDPDSAnalyzer::beginJob() } + for(size_t oc=0; oc _flash_time_to_match; + std::vector _crttrack_time_to_match; + + art::Handle< std::vector > opflashListHandle; + // Read OpFlash Time + for (size_t s = 0; s < fOpFlashesModuleLabel.size(); s++) { + e.getByLabel(fOpFlashesModuleLabel[s], opflashListHandle); + if(!opflashListHandle.isValid()){ + std::cout << "OpFlash with label " << fOpFlashesModuleLabel[s] << " not found..." << std::endl; + throw std::exception(); + } + if(fVerbosity>0) + std::cout << "Saving OpFlashes from " << fOpFlashesModuleLabel[s] << std::endl; + for (unsigned int i = 0; i < opflashListHandle->size(); ++i) { + // Get OpFlash + art::Ptr FlashPtr(opflashListHandle, i); + recob::OpFlash Flash = *FlashPtr; + _flash_time_to_match.push_back( Flash.AbsTime() ); + } + } + unsigned nOpFlash = _flash_time_to_match.size(); + + // Read CRTTracks + art::Handle> CRTTrackHandle; + e.getByLabel(fCRTTrackModuleLabel, CRTTrackHandle); + if(!CRTTrackHandle.isValid()){ + std::cout << "CRTTrack product " << fCRTTrackModuleLabel << " not found..." << std::endl; + throw std::exception(); + } + std::vector> CRTTrackVec; + art::fill_ptr_vector(CRTTrackVec, CRTTrackHandle); + const unsigned nTracks = CRTTrackVec.size(); + + for(unsigned i = 0; i < nTracks; ++i) + { + const auto track = CRTTrackVec[i]; + _crttrack_time_to_match.push_back(track->Ts0()); + } + + _opflash_matched.clear(); + _crt_tracks_matched.clear(); + _opflash_matched.resize(nOpFlash, false); + _crt_tracks_matched.resize(nTracks, false); + + for(size_t i=0; i<_flash_time_to_match.size();i++) + { + double _opflash_time = _flash_time_to_match[i]; + for(size_t j=0; j<_crttrack_time_to_match.size(); j++) + { + double _crt_track_time = _crttrack_time_to_match[j]/1000; + double time_diff = abs(_crt_track_time-_opflash_time); + if(time_diff< 1) + { + _opflash_matched[i] = true; + _crt_tracks_matched[j] = true; + break; + } + } + } + } + + + // --- Saving OpFlashes if(fSaveOpFlashes){ @@ -603,6 +721,9 @@ void opdet::SBNDPDSAnalyzer::analyze(art::Event const& e) _flash_zerr.clear(); _flash_x.clear(); _flash_xerr.clear(); + _flash_pe_co.clear(); + _flash_pe_unco.clear(); + _flash_pmt_ratio.clear(); _flash_ophit_time.clear(); _flash_ophit_risetime.clear(); _flash_ophit_starttime.clear(); @@ -614,6 +735,7 @@ void opdet::SBNDPDSAnalyzer::analyze(art::Event const& e) art::Handle< std::vector > opflashListHandle; + int flash_idx = -1 ; // Loop over all the OpFlash labels for (size_t s = 0; s < fOpFlashesModuleLabel.size(); s++) { @@ -631,9 +753,27 @@ void opdet::SBNDPDSAnalyzer::analyze(art::Event const& e) // Get OpFlash art::Ptr FlashPtr(opflashListHandle, i); recob::OpFlash Flash = *FlashPtr; - + // Skip flash if it's outside the save window. + flash_idx+=1; + if(Flash.AbsTime()*1000fOpFlashBeamWindow[1]) continue; + if(fSaveOnlyCRTPDSMatch && !_opflash_matched[flash_idx]) continue; if(fVerbosity>0) std::cout << " * " << _nopflash << " Time [ns]=" << 1000*Flash.AbsTime() << " PE=" << Flash.TotalPE() << std::endl; + + if(fSavePEFlavourPerFlash) + { + double _current_pe_co = 0; + double _current_pe_unco = 0; + std::vector> ophit_v = flashToOpHitAssns.at(i); + for (auto ophit : ophit_v) { + int opch = ophit->OpChannel(); + double channel_pe = ophit->PE(); + if(fPDSMap.pdType(opch)=="pmt_coated") _current_pe_co+=channel_pe; + else if(fPDSMap.pdType(opch)=="pmt_uncoated") _current_pe_unco+=channel_pe; + } + _flash_pe_co.push_back(_current_pe_co); + _flash_pe_unco.push_back(_current_pe_unco); + } _flash_id.push_back( _nopflash ); _flash_time.push_back( Flash.AbsTime() ); @@ -646,6 +786,7 @@ void opdet::SBNDPDSAnalyzer::analyze(art::Event const& e) _flash_xerr.push_back( Flash.XWidth() ); _flash_z.push_back( Flash.ZCenter() ); _flash_zerr.push_back( Flash.ZWidth() ); + _flash_pmt_ratio.push_back(GetPMTRatioData(Flash.PEs())); _nopflash++; if(fSaveOpHits){ @@ -679,6 +820,135 @@ void opdet::SBNDPDSAnalyzer::analyze(art::Event const& e) } + if(fSaveCRT) + { + _tr_entry_x.clear(); + _tr_entry_y.clear(); + _tr_entry_z.clear(); + _tr_exit_x.clear(); + _tr_exit_y.clear(); + _tr_exit_z.clear(); + _tr_start_x.clear(); + _tr_start_y.clear(); + _tr_start_z.clear(); + _tr_end_x.clear(); + _tr_end_y.clear(); + _tr_end_z.clear(); + _tr_dir_x.clear(); + _tr_dir_y.clear(); + _tr_dir_z.clear(); + _tr_ts0.clear(); + _tr_ets0.clear(); + _tr_ts1.clear(); + _tr_ets1.clear(); + _tr_pe.clear(); + _tr_length.clear(); + _tr_tof.clear(); + _tr_theta.clear(); + _tr_phi.clear(); + _tr_triple.clear(); + + // Read CRTTracks + art::Handle> CRTTrackHandle; + e.getByLabel(fCRTTrackModuleLabel, CRTTrackHandle); + if(!CRTTrackHandle.isValid()){ + std::cout << "CRTTrack product " << fCRTTrackModuleLabel << " not found..." << std::endl; + throw std::exception(); + } + std::vector> CRTTrackVec; + art::fill_ptr_vector(CRTTrackVec, CRTTrackHandle); + const unsigned nTracks = CRTTrackVec.size(); + + for(unsigned i = 0; i < nTracks; ++i) + { + + const auto track = CRTTrackVec[i]; + if(track->Ts0() < fCRTSaveWindow[0] || track->Ts0() > fCRTSaveWindow[1]) continue; + if(fSaveOnlyCRTPDSMatch && !_crt_tracks_matched[i]) continue; + + const int tpc_0=0; + const int tpc_1=1; + TVector3 EntryPoint_tpc0, ExitPoint_tpc0; + TVector3 EntryPoint_tpc1, ExitPoint_tpc1; + bool cross_tpc0 = CRTTrackCrossesAV(tpc_0, *track, EntryPoint_tpc0, ExitPoint_tpc0); + bool cross_tpc1 = CRTTrackCrossesAV(tpc_1, *track, EntryPoint_tpc1, ExitPoint_tpc1); + if(fSaveOnlyAVTracks) + { + if(!cross_tpc0 && !cross_tpc1) continue; + if(cross_tpc0){ + _tr_entry_x.push_back(EntryPoint_tpc0.X()); + _tr_entry_y.push_back(EntryPoint_tpc0.Y()); + _tr_entry_z.push_back(EntryPoint_tpc0.Z()); + _tr_exit_x.push_back(ExitPoint_tpc0.X()); + _tr_exit_y.push_back(ExitPoint_tpc0.Y()); + _tr_exit_z.push_back(ExitPoint_tpc0.Z()); + } + else{ + _tr_entry_x.push_back(EntryPoint_tpc1.X()); + _tr_entry_y.push_back(EntryPoint_tpc1.Y()); + _tr_entry_z.push_back(EntryPoint_tpc1.Z()); + _tr_exit_x.push_back(ExitPoint_tpc1.X()); + _tr_exit_y.push_back(ExitPoint_tpc1.Y()); + _tr_exit_z.push_back(ExitPoint_tpc1.Z()); + } + } + const geo::Point_t start = track->Start(); + _tr_start_x.push_back(start.X()); + _tr_start_y.push_back(start.Y()); + _tr_start_z.push_back(start.Z()); + + const geo::Point_t end = track->End(); + _tr_end_x.push_back(end.X()); + _tr_end_y.push_back(end.Y()); + _tr_end_z.push_back(end.Z()); + + const geo::Vector_t dir = track->Direction(); + _tr_dir_x.push_back(dir.X()); + _tr_dir_y.push_back(dir.Y()); + _tr_dir_z.push_back(dir.Z()); + + _tr_ts0.push_back(track->Ts0()); + _tr_ets0.push_back(track->Ts0Err()); + _tr_ts1.push_back(track->Ts1()); + _tr_ets1.push_back(track->Ts1Err()); + _tr_pe.push_back(track->PE()); + _tr_length.push_back(track->Length()); + _tr_tof.push_back(track->ToF()); + _tr_theta.push_back(TMath::RadToDeg() * track->Theta()); + _tr_phi.push_back(TMath::RadToDeg() * track->Phi()); + _tr_triple.push_back(track->Triple()); + } + } + + // --- Saving SPECTDC + + if(fSaveSPECTDC) + { + art::Handle> tdcHandle; + e.getByLabel(fSPECTDCLabel, tdcHandle); + if (!tdcHandle.isValid() || tdcHandle->size() == 0){ + std::cout << "No SPECTDC products found. Skip this event." << std::endl; + return; + } + else{ + const std::vector tdc_v(*tdcHandle); + for (size_t i=0; i opdet::SBNDPDSAnalyzer::GetAllHitsTruthMatch(art::Eve return allHitsTruthMap; } + +double opdet::SBNDPDSAnalyzer::GetPMTRatioData(std::vector PE_v){ + + std::map BoxMap_PECoated; + std::map BoxMap_PEUncoated; + std::map BoxMap_NCoatedCh; + std::map BoxMap_NUncoatedCh; + + for(size_t oc=0; oc=1){ + PECoated+= BoxMap_PECoated[boxID]; + PEUncoated+=BoxMap_PEUncoated[boxID]; + TotalPE += PECoated + PEUncoated; + RatioPerBox = (PEUncoated/PECoated)*BoxMap_NCoatedCh[boxID]; + RatioPerBoxWeight += RatioPerBox * (PECoated + PEUncoated); + } + } + if(TotalPE!=0) pmtratio = RatioPerBoxWeight/TotalPE; + + + return pmtratio; +} + +// -------- Function to chekc if a CRTTrack crosses the AV of TPC(0/1) and return then entry/exit point -------- +bool opdet::SBNDPDSAnalyzer::CRTTrackCrossesAV(const int tpc, const sbnd::crt::CRTTrack &crttrack, TVector3& EntryPoint, TVector3& ExitPoint){ + + TVector3 AV_min, AV_max; + // Define the AV limits for TPC0 and TPC1 + if(tpc==0){ + AV_min = {-200, -200, 0}; + AV_max = {0, 200, 500}; + } + else if(tpc==1){ + AV_min = {0, -200, 0}; + AV_max = {200, 200, 500}; + } + else std::runtime_error("Invalid TPC number. Only 0 and 1 are valid."); + + const geo::Point_t start = crttrack.Start(); + const geo::Point_t end = crttrack.End(); + + TVector3 p1 = {start.X(), start.Y(), start.Z()}; + TVector3 p2 = {end.X(), end.Y(), end.Z()}; + + TVector3 d = p2 - p1; + for (int i = 0; i < 3; ++i) { + if (d[i] == 0.0) + d[i] = 1e-10; + } + + std::array t_min, t_max, t1, t2; + for (int i = 0; i < 3; ++i) { + t_min[i] = (AV_min[i] - p1[i]) / d[i]; + t_max[i] = (AV_max[i] - p1[i]) / d[i]; + t1[i] = std::min(t_min[i], t_max[i]); + t2[i] = std::max(t_min[i], t_max[i]); + } + + double t_entry = *std::max_element(t1.begin(), t1.end()); + double t_exit = *std::min_element(t2.begin(), t2.end()); + + if (t_entry <= t_exit && t_exit >= 0.0 && t_entry <= 1.0) { + TVector3 punto_entrada = p1 + d * t_entry; + TVector3 punto_salida = p1 + d * t_exit; + //std::cout << " Entry point " << punto_entrada.X() << " " << punto_entrada.Y() << " " << punto_entrada.Z() << std::endl; + //std::cout << " Exit point " << punto_salida.X() << " " << punto_salida.Y() << " " << punto_salida.Z() << std::endl; + EntryPoint = punto_entrada; + ExitPoint = punto_salida; + return true; + } + else + { + EntryPoint={-9999., -9999., -9999.}; + ExitPoint={-9999., -9999., -9999.}; + return false; + } + +} diff --git a/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.hh b/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.hh index f6e408416..a37e993bf 100644 --- a/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.hh +++ b/sbndcode/OpDetAnalyzer/PDSAnalyzer/SBNDPDSAnalyzer_module.hh @@ -47,11 +47,14 @@ #include "lardataobj/RecoBase/SpacePoint.h" #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/PFParticle.h" +// CRT +#include "sbnobj/SBND/CRT/CRTTrack.hh" // Cosmic rejection includes #include "sbnobj/Common/Reco/OpT0FinderResult.h" #include "sbnobj/Common/Reco/SimpleFlashMatchVars.h" #include "sbnobj/Common/Reco/CRUMBSResult.h" +#include "sbnobj/SBND/Timing/DAQTimestamp.hh" #include "lardataobj/AnalysisBase/T0.h" // Geometry and mapping @@ -115,6 +118,10 @@ private: std::map GetAllHitsTruthMatch(art::Event const& e, const std::vector > &allHits); + double GetPMTRatioData(std::vector ); + + bool CRTTrackCrossesAV(const int , const sbnd::crt::CRTTrack& , TVector3& , TVector3&); + // TTree saving options bool fSaveMCTruth; bool fSaveMCParticles; @@ -123,9 +130,15 @@ private: bool fSaveRawWaveforms; bool fSaveDeconvolvedWaveforms; bool fSaveOpHits; + bool fSaveOnlyFlashHits; bool fSaveOpFlashes; + bool fSaveCRT; + bool fSaveSPECTDC; + bool fSaveOnlyCRTPDSMatch; + bool fSaveOnlyAVTracks; bool fSaveCosmicId; - + bool fSavePEFlavourPerFlash; + // Configuration parameters int fVerbosity; bool fMakePerTrackTree; @@ -145,20 +158,31 @@ private: std::string fDeconvolvedWaveformsModuleLabel; std::vector fOpHitsModuleLabel; std::vector fOpFlashesModuleLabel; + std::string fCRTTrackModuleLabel; + std::string fSPECTDCLabel; std::string fHitsLabel; std::string fReco2Label; std::string fCosmicIdModuleLabel; std::string fOpT0FinderModuleLabel; std::string fSimpleFlashMatchModuleLabel; + // Fiducial volume for MC Particles std::vector fG4BufferBoxX; std::vector fG4BufferBoxY; std::vector fG4BufferBoxZ; std::vector fG4BeamWindow; + //OpFlash save window + std::vector fOpFlashBeamWindow={-1500000, 1500000}; + + // CRT save window + std::vector fCRTSaveWindow={-1500000, 1500000}; + // PDS mapping and geometry opdet::sbndPDMapAlg fPDSMap; + std::set fPDSBoxIDs; + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); static constexpr double fDefaultSimIDE = -999.; @@ -250,6 +274,9 @@ private: std::vector _flash_x; std::vector _flash_xerr; std::vector _flash_tpc; + std::vector _flash_pe_co; + std::vector _flash_pe_unco; + std::vector _flash_pmt_ratio; std::vector> _flash_ophit_time; std::vector> _flash_ophit_risetime; std::vector> _flash_ophit_starttime; @@ -259,6 +286,44 @@ private: std::vector> _flash_ophit_pe; std::vector> _flash_ophit_ch; + + // Saving CRT information + std::vector _tr_entry_x; + std::vector _tr_entry_y; + std::vector _tr_entry_z; + std::vector _tr_exit_x; + std::vector _tr_exit_y; + std::vector _tr_exit_z; + std::vector _tr_start_x; + std::vector _tr_start_y; + std::vector _tr_start_z; + std::vector _tr_end_x; + std::vector _tr_end_y; + std::vector _tr_end_z; + std::vector _tr_dir_x; + std::vector _tr_dir_y; + std::vector _tr_dir_z; + std::vector _tr_ts0; + std::vector _tr_ets0; + std::vector _tr_ts1; + std::vector _tr_ets1; + std::vector _tr_pe; + std::vector _tr_length; + std::vector _tr_tof; + std::vector _tr_theta; + std::vector _tr_phi; + std::vector _tr_triple; + + + // Saving SPECTDC + uint64_t event_trigger_time; + uint64_t rwm_time; + + // Vectors to store CRT/PDS matched + std::vector _crt_tracks_matched; + std::vector _opflash_matched; + + // Cosmic ID std::vector _CRUMBSScore; std::vector _sliceOrigin; diff --git a/sbndcode/OpDetAnalyzer/PDSAnalyzer/job/SBNDPDSAna_config.fcl b/sbndcode/OpDetAnalyzer/PDSAnalyzer/job/SBNDPDSAna_config.fcl index 1a836f0ff..524f81b1d 100644 --- a/sbndcode/OpDetAnalyzer/PDSAnalyzer/job/SBNDPDSAna_config.fcl +++ b/sbndcode/OpDetAnalyzer/PDSAnalyzer/job/SBNDPDSAna_config.fcl @@ -13,7 +13,13 @@ sbndPDSAna: SaveDeconvolvedWaveforms: "false" SaveOpHits: "true" SaveOpFlashes: "true" + SaveSPECTDC: "true" + SaveCRT: "true" SaveCosmicId: "false" + SaveOnlyCRTPDSMatch: "false" + SavePEFlavourPerFlash: "true" + SaveOnlyFlashHits: "true" + SaveOnlyAVTracks: "true" ## Configuration parameters Verbosity: 0 @@ -34,6 +40,8 @@ sbndPDSAna: DeconvolvedWaveformsModuleLabel: "opdecopmt" OpHitsModuleLabel: ["ophitpmt"] OpFlashesModuleLabel: ["opflashtpc0", "opflashtpc1"] + CRTTrackModuleLabel: "crttracks" + SPECTDCLabel: "tdcdecoder" HitsLabel: "gaushit" Reco2Label: "pandora" CosmicIdModuleLabel: "crumbs" @@ -46,6 +54,8 @@ sbndPDSAna: G4BufferBoxZ: [-100, 600] #cm G4BeamWindow: [-10000, 12000] #ns + OpFlashBeamWindow: [-1500000, 1500000] #ns [Full flash window] + CRTSaveWindow: [-1500000, 1500000] #ns [Full flash window] } END_PROLOG diff --git a/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatioData_tool.cc b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatioData_tool.cc new file mode 100644 index 000000000..0122165c0 --- /dev/null +++ b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatioData_tool.cc @@ -0,0 +1,231 @@ +/////////////////////////////////////////////////////////////////////// +/// File: DriftEstimatorPMTRatio_tool.cc +/// +/// Base class: DriftEstimatorBase.hh +/// +/// Tool description: this tool estimates the drift coordinate +/// from the ratio between the #PE reconstructed for the +/// uncoated/coated PMTs. It requires a calibration curve +/// (speficied in the CalibrationFile fhicl parameter). +/// Once the drift has been estimated, the photon propagation +/// time is calculated using the VUV and VIS light group velocities +/// +/// Created by Fran Nicolas, June 2022 +//////////////////////////////////////////////////////////////////////// + +#include "fhiclcpp/types/Atom.h" +#include "art/Utilities/ToolMacros.h" +#include "art/Utilities/ToolConfigTable.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "larcore/Geometry/Geometry.h" + +#include +#include +#include + +#include "DriftEstimatorBase.hh" +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" + +//ROOT includes +#include "TProfile.h" +#include "TFile.h" + +namespace lightana +{ + class DriftEstimatorPMTRatioData : DriftEstimatorBase{ + + public: + + //Configuration parameters + struct Config { + + fhicl::Atom CalibrationFile { + fhicl::Name("CalibrationFile"), + fhicl::Comment("Filepath to the calibration ROOT file") + }; + + fhicl::Atom VGroupVUV { + fhicl::Name("VGroupVUV"), + fhicl::Comment("Group velocity for VUV photons") + }; + + fhicl::Atom VGroupVIS { + fhicl::Name("VGroupVIS"), + fhicl::Comment("Group velocity for VIS photons") + }; + + }; + + // Default constructor + explicit DriftEstimatorPMTRatioData(art::ToolConfigTable const& config); + + // Method giving the estimated drift coordinate + double GetDriftPosition(std::vector PE_v) override; + + // Method giving the photon propagation + double GetPropagationTime(double drift) override; + + // Method giving the photon propagation from PE vector + double PEToPropagationTime(std::vector PE_v); + + private: + double Interpolate(double val); + + // Input filepah with calibration curve + std::string fCalibrationFile; + + // Scintillation light group velocities + double fVGroupVUV; + double fVGroupVIS; + double fVGroupVUV_I; + + //Geo properties + double fDriftDistance; + double fVISLightPropTime; + double fKinkDistance; + + // Vectors with calibration variables + std::vector fPMTRatioCal; + std::vector fDriftCal; + int fNCalBins; + double fPMTRatio_MinVal; + double fPMTRatio_MaxVal; + + // PDS mapping + opdet::sbndPDMapAlg fPDSMap; + + std::set fPDSBoxIDs; + + }; + + DriftEstimatorPMTRatioData::DriftEstimatorPMTRatioData(art::ToolConfigTable const& config) + : fCalibrationFile { config().CalibrationFile() }, + fVGroupVUV { config().VGroupVUV() }, + fVGroupVIS { config().VGroupVIS() } + { + + // Open input file + std::string file_name; + cet::search_path sp("FW_SEARCH_PATH"); + if ( !sp.find_file(fCalibrationFile, file_name) ) + throw cet::exception("DriftEstimatorPMTRatioData") << "Calibration file " << + fCalibrationFile << " not found in FW_SEARCH_PATH\n"; + + TFile* input_file = TFile::Open(file_name.c_str(), "READ"); + TProfile * hProf_Calibration = (TProfile*)input_file->Get("PMTRatioCalibrationProfile"); + + //Fill calibration variables + fNCalBins = hProf_Calibration->GetNbinsX(); + for (int ix=1; ix<=fNCalBins; ix++){ + fPMTRatioCal.push_back( hProf_Calibration->GetBinCenter(ix) ); + fDriftCal.push_back( hProf_Calibration->GetBinContent(ix) ); + } + fPMTRatio_MinVal = hProf_Calibration->GetBinCenter(1); + fPMTRatio_MaxVal = hProf_Calibration->GetBinCenter(fNCalBins); + + input_file->Close(); + + auto const& tpc = art::ServiceHandle()->TPC(); + + fDriftDistance = tpc.DriftDistance(); + fVISLightPropTime = fDriftDistance/fVGroupVIS; + fKinkDistance = 0.5*fDriftDistance*(1-fVGroupVUV/fVGroupVIS); + + fVGroupVUV_I = 1./fVGroupVUV; + + for(size_t oc=0; oc PE_v){ + + std::map BoxMap_PECoated; + std::map BoxMap_PEUncoated; + std::map BoxMap_NCoatedCh; + std::map BoxMap_NUncoatedCh; + + for(size_t oc=0; oc=1){ + PECoated+= BoxMap_PECoated[boxID]; + PEUncoated+=BoxMap_PEUncoated[boxID]; + TotalPE += PECoated + PEUncoated; + RatioPerBox = (PEUncoated/PECoated)*BoxMap_NCoatedCh[boxID]/4; + RatioPerBoxWeight += 4*RatioPerBox * (PECoated + PEUncoated); + } + } + + if(TotalPE!=0){ + double pmtratio = RatioPerBoxWeight/TotalPE; + + double drift_distance; + if(pmtratio<=fPMTRatioCal[0]) + drift_distance=fDriftCal[0]; + else if(pmtratio>=fPMTRatioCal[fNCalBins-1]) + drift_distance=fDriftCal[fNCalBins-1]; + else + drift_distance=Interpolate(pmtratio); + return drift_distance; + } + else return fDriftCal[fNCalBins-1]; + } + + double DriftEstimatorPMTRatioData::GetPropagationTime(double drift){ + + // drift is here the X coordinate + // cathode: x=0 cm, PDS: x=200 cm + if(std::abs(drift) > fKinkDistance) + return (fDriftDistance-std::abs(drift)) * fVGroupVUV_I ; + else + return std::abs(drift) * fVGroupVUV_I + fVISLightPropTime; + } + + double DriftEstimatorPMTRatioData::PEToPropagationTime(std::vector PE_v){ + + double _drift = GetDriftPosition(PE_v); + + return GetPropagationTime(_drift); + } + + double DriftEstimatorPMTRatioData::Interpolate(double val){ + + size_t upix = std::upper_bound(fPMTRatioCal.begin(), fPMTRatioCal.end(), val)-fPMTRatioCal.begin(); + + double slope = ( fDriftCal[upix]-fDriftCal[upix] ) / ( fPMTRatioCal[upix]-fPMTRatioCal[upix-1] ); + + return fDriftCal[upix-1] + slope * ( val - fPMTRatioCal[upix-1] ); + } + +} + +DEFINE_ART_CLASS_TOOL(lightana::DriftEstimatorPMTRatioData) \ No newline at end of file diff --git a/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc index cf932338d..87a5eeb1a 100644 --- a/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc +++ b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc @@ -216,6 +216,7 @@ namespace lightana drift_distance=fDriftCal[fNCalBins-1]; else drift_distance=Interpolate(pmtratio); + return drift_distance; } else return fDriftCal[fNCalBins-1]; diff --git a/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc b/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc index 3d5190370..610ca4f74 100644 --- a/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc +++ b/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc @@ -31,9 +31,15 @@ namespace lightana{ //Configuration parameters struct Config { - fhicl::Atom Threshold { - fhicl::Name("Threshold"), - fhicl::Comment("Use channels abnove this thershold") + fhicl::Atom ThresholdY { + fhicl::Name("ThresholdY"), + fhicl::Comment("Use channels abnove this threshold (Y coord)") + }; + + + fhicl::Atom ThresholdZ { + fhicl::Name("ThresholdZ"), + fhicl::Comment("Use channels abnove this threshold (Z coord)") }; fhicl::Sequence PDTypes { @@ -65,11 +71,13 @@ namespace lightana{ private: void ResetVars(); - void GetCenter(std::map PEAcc, double& center, double& width); + void GetCenter(std::map PEAcc, double& center, double& width, double& threshold); // Fhicl configuration parameters std::vector fPDTypes; - float fThreshold; + double fThresholdY; + double fThresholdZ; + bool fNormalizeByPDType; unsigned int fWeightExp; @@ -98,7 +106,8 @@ namespace lightana{ FlashGeoThreshold::FlashGeoThreshold(art::ToolConfigTable const& config) : fPDTypes{ config().PDTypes() } - , fThreshold{ config().Threshold() } + , fThresholdY{ config().ThresholdY() } + , fThresholdZ{ config().ThresholdZ() } , fNormalizeByPDType{ config().NormalizeByPDType() } , fWeightExp{ config().WeightExp() } { @@ -188,8 +197,8 @@ namespace lightana{ } // Get YZ position of selected channels (above threshold) - GetCenter(fYPEAccumulator, Ycenter, Ywidth); - GetCenter(fZPEAccumulator, Zcenter, Zwidth); + GetCenter(fYPEAccumulator, Ycenter, Ywidth, fThresholdY); + GetCenter(fZPEAccumulator, Zcenter, Zwidth, fThresholdZ); } @@ -204,7 +213,7 @@ namespace lightana{ } } - void FlashGeoThreshold::GetCenter(std::map PEAcc, double& center, double& width){ + void FlashGeoThreshold::GetCenter(std::map PEAcc, double& center, double& width, double& threshold){ // set variables center = 0.; @@ -227,7 +236,7 @@ namespace lightana{ sum2+=weight*pe.first*pe.first; // For OpFlash center consider only channels above ceertain threshold - if(pe.second>fThreshold){ + if(pe.second>threshold){ center+=weight*pe.first; weightNormCenter+=weight; } diff --git a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl index f7385048f..18051a4c2 100644 --- a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl +++ b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl @@ -10,7 +10,8 @@ FlashGeoThreshold: { tool_type: "FlashGeoThreshold" PDTypes: ["pmt_coated", "pmt_uncoated"] - Threshold: 0.8 + ThresholdY: 0.75 + ThresholdZ: 0.45 NormalizeByPDType: true WeightExp: 2 } diff --git a/sbndcode/Timing/CMakeLists.txt b/sbndcode/Timing/CMakeLists.txt index e631294d5..29d2a9562 100644 --- a/sbndcode/Timing/CMakeLists.txt +++ b/sbndcode/Timing/CMakeLists.txt @@ -3,4 +3,6 @@ lardataobj::RawData) install_headers() # install_fhicl() -# install_source() \ No newline at end of file +# install_source() + +add_subdirectory(WaveformAlignment) diff --git a/sbndcode/Timing/SBNDRawTimingObj.h b/sbndcode/Timing/SBNDRawTimingObj.h index def710e18..239675bc7 100644 --- a/sbndcode/Timing/SBNDRawTimingObj.h +++ b/sbndcode/Timing/SBNDRawTimingObj.h @@ -37,7 +37,23 @@ namespace raw { uint16_t postPercent; // # 0-100, represents a percentage std::vector triggerTimeTag; // ns }; + + //TODO: remove board ID + class BoardAlignment { + // Associate one of these to every opdetwaveform in the board/digitizer, one per board + public: + BoardAlignment() {}; // constructor + BoardAlignment(unsigned int boardId, std::vector shift) : + boardId(boardId), shift(shift) {}; + + //TODO: remove boardId after validation check + unsigned int boardId; //board 0-7 + std::vector shift; //shift per waveform + + //TODO: add a flag if jittering is applied + //or a flag indicating the level of PMT timing + }; } } -#endif \ No newline at end of file +#endif diff --git a/sbndcode/Timing/WaveformAlignment/CMakeLists.txt b/sbndcode/Timing/WaveformAlignment/CMakeLists.txt new file mode 100644 index 000000000..9fd05d4f9 --- /dev/null +++ b/sbndcode/Timing/WaveformAlignment/CMakeLists.txt @@ -0,0 +1,17 @@ +set ( + MODULE_LIBRARIES + lardata::Utilities + lardataobj::RecoBase + art_root_io::tfile_support + art_root_io::TFileService_service + canvas::canvas + ROOT::Gpad + ROOT::Core + ROOT::Tree + sbnobj::SBND_Timing + sbnobj::SBND_CRT +) +cet_build_plugin(WaveformAlignment art::module SOURCE WaveformAlignment_module.cc LIBRARIES ${MODULE_LIBRARIES}) + +install_fhicl() +install_source() diff --git a/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc b/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc new file mode 100644 index 000000000..b7a28af33 --- /dev/null +++ b/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc @@ -0,0 +1,1135 @@ +//////////////////////////////////////////////////////////////////////// +// Class: WaveformAlignment +// Plugin Type: Producer +// File: WaveformAlignment_module.cc +// +// Author: Lan Nguyen (vclnguyen@ucsb.edu) +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art_root_io/TFileService.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "canvas/Persistency/Provenance/ProcessConfiguration.h" +#include "canvas/Persistency/Provenance/ProcessHistory.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "canvas/Persistency/Common/Assns.h" +#include "art/Persistency/Common/PtrMaker.h" + +#include "TCanvas.h" +#include "TTree.h" +#include "TH1D.h" +#include "TFile.h" +#include "TF1.h" +#include "TGraph.h" +#include "TMultiGraph.h" +#include "TLegend.h" +#include "TFitResult.h" +#include "TStyle.h" +#include "TSystem.h" + +#include "lardataobj/RawData/OpDetWaveform.h" +#include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "sbndcode/Timing/SBNDRawTimingObj.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + +#include +#include + +namespace sbnd { + class WaveformAlignment; +} + + +class sbnd::WaveformAlignment : public art::EDProducer { +public: + explicit WaveformAlignment(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + WaveformAlignment(WaveformAlignment const&) = delete; + WaveformAlignment(WaveformAlignment&&) = delete; + WaveformAlignment& operator=(WaveformAlignment const&) = delete; + WaveformAlignment& operator=(WaveformAlignment&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + + // Selected optional functions. + void beginJob() override; + void endJob() override; + + void ResetEventVars(); + + std::pair FitFtrig(art::Ptr wf, const int flashId, const int boardId); + double FindNearestTdc(const double tsFtrig, const std::vector tdcVec); + bool CheckShift(const double shift, const int boardId); + + void PlotFtrigCompare(const int flashId); + void PlotFtrigFit(TGraph *g, TF1 *fitf, const bool converged, const int boardId, const int flashId); + void ResetPlotVars(); + + double getPMTResponseTime(int ); + + +private: + + //---GLOBAL PARAMETERS + + enum TimingStatus { + kUndefined = -1, + kGood, //0 + kFitFailure, //1 + kOutOfBound //2 + }; + + // Plotting + std::stringstream _plotName; //raw waveform hist name + std::map tickVec; + std::map> xVec; + std::map> yVec; + std::map boardMidX; + std::map boardMidY; + + // Shifting + std::vector boardJitter[9]; + std::vector boardStatus[9]; + + // Flash counter + int nFtrigFlash; + int nPmtFlash; + + //---TREE PARAMETERS + TTree *fTree; + art::ServiceHandle tfs; + + std::vector pmtBoard{"0","1","2","3","4","5","6","7","8"}; + std::vector pmtCol{1, 2, 3, 4, 5, 6, 7, 28, 30}; + + int _run, _subrun, _event; + + // TDC stuff + std::vector _tdc_ch0; + std::vector _tdc_ch1; + std::vector _tdc_ch2; + std::vector _tdc_ch3; + std::vector _tdc_ch4; + + // PMT Timing + uint16_t _pmt_timing_type; + uint16_t _pmt_timing_ch; + + //---SERVICE + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; + + //---FHICL CONFIG PARAMETERS + + // Product label + art::InputTag fTdcDecodeLabel; + art::InputTag fTimingRefLabel; + art::InputTag fPmtDecodeLabel; + art::InputTag fPmtBoardLabel; + art::InputTag fTimingDecodeLabel; + art::InputTag fTimingBoardLabel; + art::InputTag fFtrigDecodeLabel; + art::InputTag fFtrigBoardLabel; + + // DAQ + double fWfLength; + std::vector fTdc3CaenOffset; + int fnPmtBoard; + int fnTimingBoard; + int fnChperBoard; + + // Shift configuration + int fTickLbPmt; + int fTickUbPmt; + + int fTickLbTiming; + int fTickUbTiming; + + int fPmtFitBound; + int fTimingFitBound; + + double fPmtJitterBound; + double fTimingJitterBound; + bool fCorrectCableOnly; + bool fCorrectPMTResponseTime; + + // Debug + bool fDebugTdc; + bool fDebugTimeRef; + bool fDebugFtrig; + bool fDebugPmt; + bool fDebugTiming; + + // New product labels + std::string fFtrigNewLabel; + std::string fFtrigBoardNewLabel; + + std::string fPmtNewLabel; + std::string fPmtBoardNewLabel; + std::string fPmtAlignNewLabel; + + std::string fTimingNewLabel; + std::string fTimingBoardNewLabel; + std::string fTimingAlignNewLabel; + + + std::string fPMTResponseFileName; + std::ifstream fPMTResponseFile; + + + // Plotting + bool fSaveGoodFit; + bool fSaveBadFit; + bool fSaveCompare; + std::string fSavePath; + + +}; + + +sbnd::WaveformAlignment::WaveformAlignment(fhicl::ParameterSet const& p) + : EDProducer{p} // + // More initializers here. +{ + fTdcDecodeLabel = p.get("TdcDecodeLabel", "tdcdecoder"); + fTimingRefLabel = p.get("fTimingRefLabel", "pmtdecoder"); + + fPmtDecodeLabel = p.get("PmtDecodeLabel", "pmtdecoder:PMTChannels"); + fPmtBoardLabel = p.get("PmtBoardLabel", "pmtdecoder:PMTTiming"); + + fTimingDecodeLabel = p.get("TimingDecodeLabel", "pmtdecoder:TimingChannels"); + fTimingBoardLabel = p.get("TimingBoardLabel", "pmtdecoder:TimingTiming"); + + fFtrigDecodeLabel = p.get("FtrigDecodeLabel", "pmtdecoder:FTrigChannels"); + fFtrigBoardLabel = p.get("FtrigBoardLabel", "pmtdecoder:FTrigTiming"); + + fWfLength = p.get("WfLength", 5000); + fTdc3CaenOffset = p.get>("Tdc3CaenOffset", {101, 101, 101, 101, 101, 101, 101, 101, 101}); + + fnPmtBoard = p.get("nPmtBoard", 8); + fnTimingBoard = p.get("nTimingBoard", 1); + + fnChperBoard = p.get("nChperBoard", 15); + + fTickLbPmt = p.get("TickLbPmt", 1005); + fTickUbPmt = p.get("TickUbPmt", 1035); + + fTickLbTiming = p.get("TickLbTiming", 975); + fTickUbTiming = p.get("TickUbTiming", 1005); + + fPmtFitBound = p.get("PmtFitBound", 35); + fTimingFitBound = p.get("TimingFitBound", 85); + + fPmtJitterBound = p.get("PmtJitterBound", 16); + fTimingJitterBound = p.get("TimingJitterBound", 68); + fCorrectCableOnly = p.get("CorrectCableOnly", false); + fCorrectPMTResponseTime = p.get("CorrectPMTResponseTime", true); + + fDebugTdc = p.get("DebugTdc", false); + fDebugTimeRef = p.get("DebugTimeRef", false); + fDebugFtrig = p.get("DebugFtrig", false); + fDebugPmt = p.get("DebugPmt", false); + fDebugTiming = p.get("DebugTiming", false); + + fSaveGoodFit = p.get("SaveGoodFit", false); + fSaveBadFit = p.get("SaveBadFit", false); + fSaveCompare = p.get("SaveCompare", false); + fSavePath = p.get("SavePath", ""); + + fFtrigNewLabel = p.get("FtrigNewLabel","FTrigChannels"); + fFtrigBoardNewLabel = p.get("FtrigBoardNewLabel","FTrigTiming"); + + fPmtNewLabel = p.get("PmtNewLabel","PMTChannels"); + fPmtBoardNewLabel = p.get("PmtBoardNewLabel","PMTTiming"); + fPmtAlignNewLabel = p.get("PmtAlignNewLabel","PMTAlign"); + + fTimingNewLabel = p.get("TimingNewLabel","TimingChannels"); + fTimingBoardNewLabel = p.get("TimingBoardNewLabel","TimingTiming"); + fTimingAlignNewLabel = p.get("TimingAlignNewLabel","TimingAlign"); + + fPMTResponseFileName = p.get("PMTResponseFileName",""); + + produces< std::vector< raw::OpDetWaveform > >(fFtrigNewLabel); + produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(fFtrigBoardNewLabel); + + produces< std::vector< raw::OpDetWaveform > >(fPmtNewLabel); + produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(fPmtBoardNewLabel); + produces< art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform > >(fPmtAlignNewLabel); + + produces< std::vector< raw::OpDetWaveform > >(fTimingNewLabel); + produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(fTimingBoardNewLabel); + produces< art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform > >(fTimingAlignNewLabel); + + produces< std::vector< raw::pmt::BoardAlignment > >(); +} + +void sbnd::WaveformAlignment::produce(art::Event& e) +{ + // Output data products + + // FTRIG + std::unique_ptr< std::vector< raw::OpDetWaveform > > newFtrigWf (new std::vector< raw::OpDetWaveform >); + art::PtrMaker make_ftrigwf_ptr{e, fFtrigNewLabel}; + + std::unique_ptr< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform> > newFtrigBoardAssn (new art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform >); + + // PMT + std::unique_ptr< std::vector< raw::OpDetWaveform > > newPmtWf (new std::vector< raw::OpDetWaveform >); + art::PtrMaker make_pmtwf_ptr{e, fPmtNewLabel}; + + std::unique_ptr< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform> > newPmtBoardAssn (new art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform >); + std::unique_ptr< art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform> > newPmtAlignAssn (new art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform >); + + // TIMING + std::unique_ptr< std::vector< raw::OpDetWaveform > > newTimingWf (new std::vector< raw::OpDetWaveform >); + art::PtrMaker make_timingwf_ptr{e, fTimingNewLabel}; + + std::unique_ptr< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform> > newTimingBoardAssn (new art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform >); + std::unique_ptr< art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform> > newTimingAlignAssn (new art::Assns< raw::pmt::BoardAlignment, raw::OpDetWaveform >); + + // Board Alignment + std::unique_ptr< std::vector< raw::pmt::BoardAlignment >> newBoardAlign (new std::vector< raw::pmt::BoardAlignment >); + art::PtrMaker make_align_ptr{e}; + + ResetEventVars(); + + if(fSaveGoodFit | fSaveBadFit | fSaveCompare){ + gSystem->Exec(Form("mkdir -p %s", fSavePath.c_str())); + + gStyle->SetFrameLineWidth(2); + gStyle->SetTextFont(62); + gStyle->SetTextSize(0.07); + gStyle->SetLabelFont(62, "xyz"); + gStyle->SetLabelSize(0.04, "xyz"); + gStyle->SetTitleSize(0.04, "xyz"); + gStyle->SetTitleFont(62, "xyz"); + gStyle->SetTitleX(0.55); //title X location + gStyle->SetTitleY(1.1); //title Y location + gStyle->SetTitleW(0.75); //title width + gStyle->SetTitleH(0.3); //title height + gStyle->SetTitleFontSize(0.04); + } + + _run = e.id().run(); + _subrun = e.id().subRun(); + _event = e.id().event(); + + if (fDebugTdc | fDebugTimeRef | fDebugFtrig | fDebugPmt | fDebugTiming) + std::cout <<"#----------RUN " << _run << " SUBRUN " << _subrun << " EVENT " << _event <<"----------#\n"; + + //------------------------PMT Timing--------------------------// + art::Handle timingRefHandle; + e.getByLabel(fTimingRefLabel, timingRefHandle); + + if (!timingRefHandle.isValid()){ + throw cet::exception("WaveformAlignment") << "No raw::TimingReferenceInfo found w/ tag " << fTimingRefLabel << ". Check data quality!"; + } + else{ + raw::TimingReferenceInfo const& pmt_timing(*timingRefHandle); + + _pmt_timing_type = pmt_timing.timingType; + _pmt_timing_ch = pmt_timing.timingChannel; + + if (fDebugTimeRef){ + std::cout << "Timing Reference For Decoding PMT" << std::endl; + std::cout << " Type = " << _pmt_timing_type << " (SPECTDC = 0; PTB HLT = 1; CAEN-only = 3)." << std::endl; + std::cout << " Channel = " << _pmt_timing_ch << " (TDC ETRIG = 4; PTB BNB Beam+Light = 2)." << std::endl; + } + } + timingRefHandle.removeProduct(); + + //---------------------------TDC-----------------------------// + art::Handle> tdcHandle; + e.getByLabel(fTdcDecodeLabel, tdcHandle); + + if (!tdcHandle.isValid() || tdcHandle->size() == 0){ + throw cet::exception("WaveformAlignment") << "No sbnd::timing::DAQTimestamp found w/ tag " << fTdcDecodeLabel << ". Check data quality!"; + } + else{ + + std::vector tdc_v(*tdcHandle); + + for (size_t i=0; i ftrigBoardAssn(wf_ftrig_v, e, fFtrigBoardLabel); + + int nTotalBoard = fnPmtBoard + fnTimingBoard; + + if (std::fmod(wf_ftrig_v.size(), nTotalBoard) != 0) + throw cet::exception("WaveformAlignment") << "raw::OpDetWaveform found w/ tag " << fFtrigDecodeLabel << " does not have the same number of flashes per board. Check data quality!"; + + nFtrigFlash = wf_ftrig_v.size() / nTotalBoard; + + if (fDebugFtrig) + std::cout << std::endl << "Found OpDetWaveform FTRIG size = " << wf_ftrig_v.size() << ", nFTRIG per board = " << nFtrigFlash << std::endl; + + for (int flashIdx = 0; flashIdx < nFtrigFlash; flashIdx++){ + + if (fDebugFtrig) std::cout << " Looping over OpDet id " << flashIdx << "..." << std::endl; + + for (int boardIdx = 0; boardIdx < nTotalBoard; boardIdx++){ //wf->ChannelNumber() and j are the same + + int status = kUndefined; + + int wfIdx = flashIdx + boardIdx * nFtrigFlash; + art::Ptr wf(wf_ftrig_v.at(wfIdx)); + + std::vector> wf_board_v(ftrigBoardAssn.at(wf.key())); + if(wf_board_v.size() != 1 ) + throw cet::exception("WaveformAlignment") << "No raw::wf::BoardTimingInfo found w/ tag" << fFtrigBoardLabel <<". Check data quality!"; + + art::Ptr wf_board(wf_board_v.front()); + + //Fit FTRIG rising edge + std::pair midPoint = FitFtrig(wf, flashIdx, wf->ChannelNumber()); + double tickFtrig = midPoint.first; + double shift = 0; //Default as 0 + + //Check if fit converged + if (tickFtrig == std::numeric_limits::min()){ + status = kFitFailure; + if (fDebugFtrig) std::cout << " board id = " << wf->ChannelNumber() << " has FTRIG fit failure." << std::endl; + }else{ + + //convert tick to ts + double tsFtrig = wf_board->triggerTimeTag[0] - (fWfLength - tickFtrig)*2.0; + + //TDC as reference frame TODO: add picosecond correction + if ((_pmt_timing_type == 0) & (_pmt_timing_ch == 4)){ + + double tsFrame = FindNearestTdc(tsFtrig, _tdc_ch3); + + double tempShift = fTdc3CaenOffset[wf->ChannelNumber()] - std::abs(tsFtrig - tsFrame); + + if (CheckShift(tempShift, wf->ChannelNumber())){ + status = kGood; + + shift = std::round(tempShift * 1000.0) / 1000.0; //round to 3 decimal place of ns + + //shift = tempShift; + if (fDebugFtrig) std::cout << " board id = " << wf->ChannelNumber() << " has rising tick value = " << tickFtrig << " and shift value = " << shift << " ns." << std::endl; + }else{ + status = kOutOfBound; //either jitter too much or no reference frame + } + } + + //TODO: Add reference for PTB and CAEN only + + } + + //Save to vector + boardJitter[wf->ChannelNumber()].push_back(shift); + boardStatus[wf->ChannelNumber()].push_back(status); + + //Save a fraction of FTRIG waveform + int saveLb = std::numeric_limits::min(); + int saveUb = std::numeric_limits::min(); + if (wf->ChannelNumber() == 8){ + saveLb = fTickLbTiming - 25; + saveUb = fTickLbTiming + 100; + }else{ + saveLb = fTickLbPmt - 25; + saveUb = fTickLbPmt + 100; + } + double saveTs = wf->TimeStamp() + saveLb*2.0; + std::vector adc_vec(saveUb - saveLb, 0); + + for (int i = 0; i < (saveUb - saveLb); i++){ + adc_vec[i] = wf->Waveform()[i+saveLb]; + } + + raw::OpDetWaveform new_wf(saveTs, wf->ChannelNumber(), adc_vec); + newFtrigWf->push_back(new_wf); + art::Ptr wfPtr = make_ftrigwf_ptr(newFtrigWf->size()-1); + newFtrigBoardAssn->addSingle(wf_board, wfPtr); + + if (fSaveCompare){ + //boardMidX[wf->ChannelNumber()] = tsFtrig; + //boardMidY[wf->ChannelNumber()] = midPoint.second; + tickVec[wf->ChannelNumber()] = tickFtrig; + for (size_t i = 0; i < wf->Waveform().size(); i++){ + yVec[wf->ChannelNumber()].push_back((double)wf->Waveform()[i]); + xVec[wf->ChannelNumber()].push_back(wf_board->triggerTimeTag[0] - (fWfLength - i)*2.0); + } + } + } //Done looping over boards + + if (fSaveCompare){ + PlotFtrigCompare(flashIdx); + ResetPlotVars(); + } + } //Done looping over OpDetWf + + // Save board alignment product + for (size_t i = 0; i < std::size(boardJitter); i++){ + raw::pmt::BoardAlignment board_align; + board_align.boardId = i; + + std::vector shift_vec; + for (auto j: boardJitter[i]){ + shift_vec.push_back(j); + } + board_align.shift = shift_vec; + newBoardAlign->push_back(board_align); + } + } + wfFtrigHandle.removeProduct(); + wf_ftrig_v.clear(); + + //---------------------------ALIGN PMT WAVEFORMS-----------------------------// + art::Handle> wfPmtHandle; + std::vector> wf_pmt_v; + e.getByLabel(fPmtDecodeLabel, wfPmtHandle); + + if (!wfPmtHandle.isValid() || wfPmtHandle->size() == 0){ + throw cet::exception("WaveformAlignment") << "No raw::OpDetWaveform found w/ tag " << fPmtDecodeLabel << ". Check data quality!"; + } + else{ + + art::fill_ptr_vector(wf_pmt_v, wfPmtHandle); + art::FindManyP pmtBoardAssn(wf_pmt_v, e, fPmtBoardLabel); + + if (std::fmod((wf_pmt_v.size() / fnChperBoard), fnPmtBoard) != 0) + throw cet::exception("WaveformAlignment") << "raw::OpDetWaveform found w/ tag " << fPmtDecodeLabel << " does not have the same number of flashes per board. Check data quality!"; + nPmtFlash = (wf_pmt_v.size() / fnChperBoard) / fnPmtBoard; + + if (fDebugPmt) std::cout << std::endl << "Found OpDetWaveform PMT size = " << wf_pmt_v.size() << ", nPmtFlash per board = " << nPmtFlash << std::endl; + + for (int boardIdx = 0 ; boardIdx < fnPmtBoard; boardIdx++){ + + if (fDebugPmt) std::cout << " Looping over board " << boardIdx << "..." << std::endl; + + for (int flashIdx = 0; flashIdx < nPmtFlash; flashIdx++){ + + if (fDebugPmt) std::cout << " flash " << flashIdx << " has shift value " << boardJitter[boardIdx][flashIdx] << std::endl; + for (int chIdx = 0; chIdx < fnChperBoard; chIdx++){ + + int wfIdx = boardIdx * nPmtFlash * fnChperBoard + flashIdx * fnChperBoard + chIdx; + art::Ptr wf(wf_pmt_v.at(wfIdx)); + + //Get assn + std::vector> wf_board_v(pmtBoardAssn.at(wf.key())); + if(wf_board_v.size() != 1 ) + throw cet::exception("WaveformAlignment") << "No raw::wf::BoardTimingInfo found w/ tag" << fPmtBoardLabel <<". Check data quality!"; + + art::Ptr wf_board(wf_board_v.front()); + art::Ptr wf_align = make_align_ptr(boardIdx); + + fPMTCalibrationDatabaseService = lar::providerFrom(); + double total_transit = fPMTCalibrationDatabaseService->getTotalTransitTime(wf->ChannelNumber()); + + double correction = 0; + correction -= total_transit; //total transit = propagation time of photon from PMT to digitiser + if(fCorrectCableOnly == false) correction -= boardJitter[boardIdx][flashIdx]; + if(fCorrectPMTResponseTime) + correction-=getPMTResponseTime(wf->ChannelNumber()); + correction /= 1000; //ns to us conversion + + double new_ts = wf->TimeStamp() + correction; //us + //std::cout << std::setprecision(6) << correction << std::endl; + //std::cout << std::setprecision(9)<< "old ts = " << wf->TimeStamp() << ", new ts = " << new_ts << std::endl; + + std::vector adc_vec(wf->Waveform().size(), 0); + for (size_t i = 0; i < wf->Waveform().size(); i++){ + adc_vec[i] = wf->Waveform()[i]; + } + + raw::OpDetWaveform new_wf(new_ts, wf->ChannelNumber(), adc_vec); + + newPmtWf->push_back(new_wf); + art::Ptr wfPtr = make_pmtwf_ptr(newPmtWf->size()-1); + + newPmtBoardAssn->addSingle(wf_board, wfPtr); + newPmtAlignAssn->addSingle(wf_align, wfPtr); + } + } + } + } + wfPmtHandle.removeProduct(); + wf_pmt_v.clear(); + + //---------------------------ALIGN TIMING WAVEFORMS-----------------------------// + art::Handle> wfTimingHandle; + std::vector> wf_timing_v; + e.getByLabel(fTimingDecodeLabel, wfTimingHandle); + + if (!wfTimingHandle.isValid() || wfTimingHandle->size() == 0){ + throw cet::exception("WaveformAlignment") << "No raw::OpDetWaveform found w/ tag " << fTimingDecodeLabel << ". Check data quality!"; + } + else{ + + art::fill_ptr_vector(wf_timing_v, wfTimingHandle); + art::FindManyP timingBoardAssn(wf_timing_v, e, fTimingBoardLabel); + + //There is only 1 timing CAEN in the hardware and not every channel is saved + int nChTiming = wf_timing_v.size()/nPmtFlash; + + if (fDebugTiming) std::cout << std::endl << "Found OpDetWaveform Timing size = " << wf_timing_v.size() << ", nFlash per board = " << nPmtFlash << std::endl; + + //TODO: fix hardcoded + int boardIdx = 8; + + if (fDebugTiming) std::cout << " Looping over board " << boardIdx << "..." << std::endl; + + for (int flashIdx = 0; flashIdx < nPmtFlash; flashIdx++){ + + if (fDebugTiming) std::cout << " flash " << flashIdx << " has shift value " << boardJitter[boardIdx][flashIdx] << std::endl; + + for (int chIdx = 0; chIdx < nChTiming; chIdx++){ + + int wfIdx = flashIdx * nChTiming + chIdx; + art::Ptr wf(wf_timing_v.at(wfIdx)); + + //Get assn + std::vector> wf_board_v(timingBoardAssn.at(wf.key())); + if(wf_board_v.size() != 1 ) + throw cet::exception("WaveformAlignment") << "No raw::wf::BoardTimingInfo found w/ tag" << fTimingBoardLabel <<". Check data quality!"; + + art::Ptr wf_board(wf_board_v.front()); + art::Ptr wf_align = make_align_ptr(boardIdx); + + //Timing CAEN is not connected to PMT + double correction = boardJitter[boardIdx][flashIdx]; + correction /= 1000; //ns to us + + double new_ts = wf->TimeStamp() + correction; //ns to us + std::vector adc_vec(wf->Waveform().size(), 0); + + for (size_t i = 0; i < wf->Waveform().size(); i++){ + adc_vec[i] = wf->Waveform()[i]; + } + + raw::OpDetWaveform new_wf(new_ts, wf->ChannelNumber(), adc_vec); + + newTimingWf->push_back(new_wf); + art::Ptr wfPtr = make_timingwf_ptr(newTimingWf->size()-1); + + newTimingBoardAssn->addSingle(wf_board, wfPtr); + newTimingAlignAssn->addSingle(wf_align, wfPtr); + } + } + } + wfTimingHandle.removeProduct(); + wf_timing_v.clear(); + + if (fDebugTdc | fDebugTimeRef | fDebugFtrig | fDebugPmt | fDebugTiming) + std::cout <<"#--------------------------------------------------------#" << std::endl; + + //Put product in event + e.put(std::move(newBoardAlign)); + + e.put(std::move(newFtrigWf), fFtrigNewLabel); + e.put(std::move(newFtrigBoardAssn), fFtrigBoardNewLabel); + + e.put(std::move(newPmtWf), fPmtNewLabel); + e.put(std::move(newPmtBoardAssn), fPmtBoardNewLabel); + e.put(std::move(newPmtAlignAssn), fPmtAlignNewLabel); + + e.put(std::move(newTimingWf), fTimingNewLabel); + e.put(std::move(newTimingBoardAssn), fTimingBoardNewLabel); + e.put(std::move(newTimingAlignAssn), fTimingAlignNewLabel); + + //Fill once every event + fTree->Fill(); +} + +void sbnd::WaveformAlignment::beginJob() +{ + //Event Tree + fTree = tfs->make("events", ""); + + fTree->Branch("run", &_run); + fTree->Branch("subrun", &_subrun); + fTree->Branch("event", &_event); + fTree->Branch("tdc_ch0", &_tdc_ch0); + fTree->Branch("tdc_ch1", &_tdc_ch1); + fTree->Branch("tdc_ch2", &_tdc_ch2); + fTree->Branch("tdc_ch3", &_tdc_ch3); + fTree->Branch("tdc_ch4", &_tdc_ch4); + + fTree->Branch("pmt_timing_type", &_pmt_timing_type); + fTree->Branch("pmt_timing_ch", &_pmt_timing_ch); + + fTree->Branch("shift_board_0", &boardJitter[0]); + fTree->Branch("shift_board_1", &boardJitter[1]); + fTree->Branch("shift_board_2", &boardJitter[2]); + fTree->Branch("shift_board_3", &boardJitter[3]); + fTree->Branch("shift_board_4", &boardJitter[4]); + fTree->Branch("shift_board_5", &boardJitter[5]); + fTree->Branch("shift_board_6", &boardJitter[6]); + fTree->Branch("shift_board_7", &boardJitter[7]); + fTree->Branch("shift_board_8", &boardJitter[8]); + + fTree->Branch("status_board_0", &boardStatus[0]); + fTree->Branch("status_board_1", &boardStatus[1]); + fTree->Branch("status_board_2", &boardStatus[2]); + fTree->Branch("status_board_3", &boardStatus[3]); + fTree->Branch("status_board_4", &boardStatus[4]); + fTree->Branch("status_board_5", &boardStatus[5]); + fTree->Branch("status_board_6", &boardStatus[6]); + fTree->Branch("status_board_7", &boardStatus[7]); + fTree->Branch("status_board_8", &boardStatus[8]); +} + +void sbnd::WaveformAlignment::endJob() +{ +} + +void sbnd::WaveformAlignment::ResetEventVars() +{ + _run = -1; _subrun = -1; _event = -1; + + _tdc_ch0.clear(); + _tdc_ch1.clear(); + _tdc_ch2.clear(); + _tdc_ch3.clear(); + _tdc_ch4.clear(); + + _pmt_timing_type = -1; + _pmt_timing_ch = -1; + + nFtrigFlash = -1; + nPmtFlash = -1; + + for (size_t i = 0; i < std::size(boardJitter); ++i){ + boardJitter[i].clear(); + boardStatus[i].clear(); + } +} + +std::pair sbnd::WaveformAlignment::FitFtrig(art::Ptr wf, const int flashId, const int boardId) +{ + //--------Make TGraph + size_t wfLength = wf->Waveform().size(); + std::vector x(wfLength), y(wfLength); + + for (size_t i = 0; i < wfLength; i++){ + y[i] = (double)wf->Waveform()[i]; + x[i] = i; + } + + TGraph *g = new TGraph(wfLength, &x[0], &y[0]); + + //---------Find baseline value + double baseline = 0.; + for(unsigned int i = 0 ; i < 800; i++) { + baseline += y[i]; + } + baseline /= 800; + + //-------Guess rising tick + double tickLb = std::numeric_limits::min(); + double tickUb = std::numeric_limits::min(); + + if (boardId == 8){ + tickLb = fTickLbTiming; + tickUb = fTickUbTiming; + }else{ + tickLb = fTickLbPmt; + tickUb = fTickUbPmt; + } + + double largestRise = std::numeric_limits::min(); + int risingTickGuess = std::numeric_limits::min(); + for(int i = tickLb; i < tickUb; i++){ + + double diff = y[i + 1] - y[i]; + if(diff > largestRise){ + largestRise = diff; + risingTickGuess = i; + } + } + + //-------Find amplitude + double amplitude = 0; + for(int i = risingTickGuess; i < risingTickGuess+50; i++){ + amplitude += y[i]; + } + amplitude /= 50; + amplitude -= baseline; + + //---------Define Rising Edge + + int fitLb = std::numeric_limits::min(); + int fitUb = std::numeric_limits::min(); + + if (boardId == 8){ + //TODO: Maybe need to change the equation for Timing CAEN? + fitLb = risingTickGuess - 85; + fitUb = risingTickGuess + 55; + }else{ + fitLb = risingTickGuess - fPmtFitBound; //25; + fitUb = risingTickGuess + fPmtFitBound; //25; + } + + g->GetXaxis()->SetRangeUser(x[fitLb], x[fitUb]); + + std::string rising_edge = "1/(1+exp(([0]-x)/[1]))*[2]+[3]"; + + TF1 *fitf = new TF1("fitf", rising_edge.c_str(), x[fitLb], x[fitUb]); + fitf->SetParameter(0, x[risingTickGuess]); + fitf->SetParameter(1, 0.3); + fitf->SetParameter(2, amplitude); + fitf->SetParameter(3, baseline); + + //---------Fit + TFitResultPtr fp = g->Fit(fitf,"SRQ","", x[fitLb], x[fitUb]); + bool converged = !(bool)(int(fp)); + + //-------Check if fit converges + std::pair midPoint(std::numeric_limits::min(), std::numeric_limits::min()); + if(converged) { + double parA = fitf->GetParameter(0); + double parMag = fitf->GetParameter(2); + double parOffset = fitf->GetParameter(3); + + double midPointY = parMag/2 + parOffset; + double midPointX = fitf->GetX(midPointY, parA - 20, parA + 20); + + midPoint.first = midPointX; + midPoint.second = midPointY; + } + + //-------Save fits to plots + if(converged && fSaveGoodFit) PlotFtrigFit(g, fitf, converged, boardId, flashId); + if(!converged && fSaveBadFit) PlotFtrigFit(g, fitf, converged, boardId, flashId); + + return midPoint; +} + +double sbnd::WaveformAlignment::FindNearestTdc(const double tsFtrig, const std::vector tdcVec) +{ + double smallestDiff = std::numeric_limits::max(); + double nearestTdc = std::numeric_limits::min(); + + for (size_t i = 0; i < tdcVec.size(); i++){ + double tdc = (double)tdcVec[i]; + double diff = std::abs(tsFtrig - tdc); + + if (diff < smallestDiff){ + smallestDiff = diff; + nearestTdc = tdc; + } + } + + return nearestTdc; +} + +bool sbnd::WaveformAlignment::CheckShift(const double shift, const int boardId) +{ + if (boardId == 8){ + if (std::abs(shift) > fTimingJitterBound){ + if (fDebugFtrig) std::cout << " board id = " << boardId << " has jittering " << shift << " > " << fTimingJitterBound << std::endl; + return false; + } + }else{ + if (std::abs(shift) > fPmtJitterBound){ + if (fDebugFtrig) std::cout << " board id = " << boardId << " has jittering " << shift << " > " << fPmtJitterBound << std::endl; + return false; + } + } + + //TODO: add more check? + + return true; +} + +double sbnd::WaveformAlignment::getPMTResponseTime(int channelNumber) +{ + + std::string fname; + cet::search_path sp("FW_SEARCH_PATH"); + sp.find_file(fPMTResponseFileName, fname); + std::ifstream fPMTResponseFile(fname.c_str(), std::ios::in); + if (!fPMTResponseFile.is_open()) { + std::cerr << "Could not open the file.\n"; + return 0; + } + + std::string line; + + while (std::getline(fPMTResponseFile, line)) { + std::stringstream linestream(line); + int channel_number; + double channel_delay; + linestream + >> channel_number + >> channel_delay; + if (channel_number == channelNumber) { + return channel_delay; + } + } + return 0.; +} + +void sbnd::WaveformAlignment::PlotFtrigFit(TGraph *g, TF1 *fitf, const bool converged, const int boardId, const int flashId) +{ + + _plotName.str(std::string()); + _plotName << "run" << _run << "_subrun" << _subrun << "_event" << _event << "_board" << boardId << "_flash" << flashId; + + int fitCol = 8; + + if (!converged){ + fitCol = 2; + _plotName << "_BADFIT"; + } + + gStyle->SetStatX(0.95); + gStyle->SetStatY(0.55); + gStyle->SetOptStat(0); + gStyle->SetOptFit(1); + gStyle->SetPadTickX(1); + gStyle->SetPadTickY(1); + + TCanvas *c = new TCanvas(_plotName.str().c_str(), _plotName.str().c_str(), 600, 500); + c->cd(); + + g->SetTitle(_plotName.str().c_str()); + g->GetXaxis()->SetTitle("Tick"); + g->GetYaxis()->SetTitle("ADC"); + + g->SetMarkerStyle(5); + g->SetMarkerSize(2); + g->SetMarkerColor(kBlack); + g->SetLineColor(kBlack); + g->SetLineWidth(1); + + g->Draw("ALP"); + + fitf->SetLineColor(fitCol); + fitf->Draw("SAME"); + + if (converged){ + TGraph *p = new TGraph(); + double halfPointY = fitf->GetParameter(2)/2 + fitf->GetParameter(3); + double risingEdge = fitf->GetX(halfPointY, fitf->GetParameter(0) - 20, fitf->GetParameter(0) + 20); + + p->SetPoint(0, risingEdge, halfPointY); + p->SetMarkerStyle(5); + p->SetMarkerSize(4); + p->SetMarkerColor(kRed); + p->Draw("Psame"); + } + + c->SetLeftMargin(0.15); + c->SetRightMargin(0.1); + c->SetBottomMargin(0.25); + + g->SetMinimum(1800); g->SetMaximum(6500); + g->GetYaxis()->SetNdivisions(505, true); + g->GetXaxis()->SetNdivisions(515, true); + g->GetXaxis()->SetTitleOffset(2.1); + g->GetXaxis()->CenterTitle(); + g->GetXaxis()->SetNoExponent(1); + + for (int i=1; i <= g->GetXaxis()->GetNlabels() + 2; i++){ + g->GetXaxis()->ChangeLabel(i, 60, 0.05, 31, -1, -1); + } + + c->SaveAs(Form("%s/%s.png", fSavePath.c_str(), _plotName.str().c_str())); +} + +void sbnd::WaveformAlignment::PlotFtrigCompare(const int flashId) +{ + auto *mg1 = new TMultiGraph(); + auto *mg2 = new TMultiGraph(); + auto *lg1 = new TLegend(0.2,0.5,0.35,0.9); + auto *lg2 = new TLegend(0.2,0.5,0.35,0.9); + int x_lb = std::abs(tickVec[0] - 2); + int x_ub = std::abs(tickVec[0] + 4); + + //---------Find baseline value + double baseline = 0.; + for(unsigned int i = 0 ; i < 800; i++) { + baseline += yVec[0][i]; + } + baseline /= 800; + + for(int k = 0; k < (fnPmtBoard + fnTimingBoard); k++){ + + std::vector x1 = xVec[k]; + std::vector y = yVec[k]; + + if ((x1.size() == 0) | (y.size() == 0)) continue; + + //-------Apply shift + double shift = boardJitter[k][flashId]; + std::vector x2; + for (size_t i = 0; i < x1.size(); i++){ + x2.push_back(x1[i] - shift); + } + + //------Shift baseline + double current_baseline = 0; + for(size_t i = 0 ; i < 800; i++) { + current_baseline += y[i]; + } + current_baseline /= 800; + + double baseline_shift = baseline - current_baseline; + for (size_t i = 0; i < y.size(); i++){ + y[i] += baseline_shift; + } + + auto *g1 = new TGraph(fWfLength, &x1[0], &y[0]); + g1->SetTitle(Form("Board_%s", pmtBoard[k].c_str())); + g1->SetMarkerStyle(5); + g1->SetMarkerSize(1); + g1->SetMarkerColor(pmtCol[k]); + g1->SetLineColor(pmtCol[k]); + + auto *g2 = new TGraph(fWfLength, &x2[0], &y[0]); + g2->SetTitle(Form("Board_%s", pmtBoard[k].c_str())); + g2->SetMarkerStyle(5); + g2->SetMarkerSize(1); + g2->SetMarkerColor(pmtCol[k]); + g2->SetLineColor(pmtCol[k]); + + //auto *p = new TGraph(); + //p->SetPoint(0, boardMidX[k] + shift, boardMidY[k] + baseline_shift); + //p->SetMarkerStyle(50); + //p->SetMarkerSize(2); + //p->SetMarkerColor(pmtCol[k]); + //mg2->Add(p); + + mg1->Add(g1); + mg2->Add(g2); + + lg1->AddEntry(g1, g1->GetTitle(), "lp"); + lg2->AddEntry(g2, g2->GetTitle(), "lp"); + } + _plotName.str(std::string()); + _plotName << "run" << _run << "_subrun" << _subrun << "_event" << _event << "_flash" << flashId; + + TCanvas *c = new TCanvas(_plotName.str().c_str(), _plotName.str().c_str(), 600, 1000); + c->Divide(1, 2, 0, 0); + + c->cd(1); + gPad->SetLeftMargin(0.2); + gPad->SetBottomMargin(0.3); + + mg1->Draw("ALP"); + + mg1->SetTitle(_plotName.str().c_str()); + mg1->GetXaxis()->SetTitle("Time [ns]"); + mg1->GetYaxis()->SetTitle("ADC"); + + mg1->SetMinimum(1800); mg1->SetMaximum(7000); + mg1->GetYaxis()->SetNdivisions(505, true); + + //mg1->GetXaxis()->SetLimits(xVec[0][x_lb+68], xVec[0][x_ub+68]); + mg1->GetXaxis()->SetLimits(xVec[0][x_lb], xVec[0][x_ub]); + mg1->GetXaxis()->SetNdivisions(310, true); + mg1->GetXaxis()->SetLabelOffset(-0.01); + mg1->GetXaxis()->SetTitleOffset(3.1); + mg1->GetXaxis()->SetNoExponent(1); + + for (int i=1; i <= mg1->GetXaxis()->GetNlabels() + 2; i++){ + mg1->GetXaxis()->ChangeLabel(i, 60, 0.04, 31, -1, -1); + } + + lg1->SetTextSize(0.03); + lg1->Draw(); + + c->Modified(); + + c->cd(2); + gPad->SetLeftMargin(0.2); + gPad->SetBottomMargin(0.3); + _plotName.str(std::string()); + _plotName << "run" << _run << "_subrun" << _subrun << "_event" << _event << "_flash" << flashId << "_aligned"; + + mg2->Draw("ALP"); + + mg2->SetTitle(_plotName.str().c_str()); + mg2->GetXaxis()->SetTitle("Time [ns]"); + mg2->GetYaxis()->SetTitle("ADC"); + + mg2->SetMinimum(1800); mg2->SetMaximum(7000); + mg2->GetYaxis()->SetNdivisions(505, true); + + //mg2->GetXaxis()->SetLimits(xVec[0][x_lb+68], xVec[0][x_ub+68]); + mg2->GetXaxis()->SetLimits(xVec[0][x_lb], xVec[0][x_ub]); + mg2->GetXaxis()->SetNdivisions(310, true); + mg2->GetXaxis()->SetLabelOffset(-0.01); + mg2->GetXaxis()->SetTitleOffset(3.1); + mg2->GetXaxis()->SetNoExponent(1); + + for (int i=1; i <= mg2->GetXaxis()->GetNlabels() + 2; i++){ + mg2->GetXaxis()->ChangeLabel(i, 60, 0.04, 31, -1, -1); + } + + lg2->SetTextSize(0.03); + lg2->Draw(); + + _plotName.str(std::string()); + _plotName << "run" << _run << "_subrun" << _subrun << "_event" << _event << "_flash" << flashId; + + c->Modified(); + + c->SaveAs(Form("%s/%s.png", fSavePath.c_str(), _plotName.str().c_str())); +} + +void sbnd::WaveformAlignment::ResetPlotVars() +{ + tickVec.clear(); + xVec.clear(); + yVec.clear(); + boardMidX.clear(); + boardMidY.clear(); +} +DEFINE_ART_MODULE(sbnd::WaveformAlignment) diff --git a/sbndcode/Timing/WaveformAlignment/run_wf_align.fcl b/sbndcode/Timing/WaveformAlignment/run_wf_align.fcl new file mode 100644 index 000000000..77df560c1 --- /dev/null +++ b/sbndcode/Timing/WaveformAlignment/run_wf_align.fcl @@ -0,0 +1,44 @@ +#include "services_sbnd.fcl" +#include "wfAlignConfig.fcl" + +process_name: WFALIGN # The process name must NOT contain any underscores + +source: +{ + module_type: RootInput # Telling art we want a ROOT input + maxEvents: -1 +} + +services: +{ + TFileService: { fileName: "wfAlignOutput.root" } + @table::sbnd_services +} + +outputs: +{ + out1: + { + fileName: "%ifb_wfalign.root" + module_type: RootOutput + dataTier: "reco" + #outputCommands: [ + # "keep *_*_*_*" + # , "drop *_pmtdecoder_*_*" + # ] + } +} + +physics: +{ + producers: + { + wfalign: @local::wfAlign + } + + prod: [ wfalign ] + stream1: [ out1 ] + + trigger_paths: [ prod ] + end_paths: [ stream1 ] +} diff --git a/sbndcode/Timing/WaveformAlignment/wfAlignConfig.fcl b/sbndcode/Timing/WaveformAlignment/wfAlignConfig.fcl new file mode 100644 index 000000000..c80c944ab --- /dev/null +++ b/sbndcode/Timing/WaveformAlignment/wfAlignConfig.fcl @@ -0,0 +1,9 @@ +BEGIN_PROLOG + +wfAlign: +{ + module_type: "WaveformAlignment" + PMTResponseFileName: "PMTResponseTime.txt" + CorrectPMTResponseTime: true +} +END_PROLOG diff --git a/sbndcode/Timing/classes_def.xml b/sbndcode/Timing/classes_def.xml index 4e7fa4847..c916b62ee 100644 --- a/sbndcode/Timing/classes_def.xml +++ b/sbndcode/Timing/classes_def.xml @@ -23,4 +23,16 @@ + + + + + + + + + + + +