diff --git a/duneana/SolarNuAna/SolarNuAna_module.cc b/duneana/SolarNuAna/SolarNuAna_module.cc index 85e07ec1..73b07807 100644 --- a/duneana/SolarNuAna/SolarNuAna_module.cc +++ b/duneana/SolarNuAna/SolarNuAna_module.cc @@ -79,7 +79,7 @@ namespace solar void ResetVariables(); // --- Our fcl parameter labels for the modules that made the data products - std::string fHitLabel, fTrackLabel, fOpHitLabel, fOpFlashLabel, fGEANTLabel; + std::string fHitLabel, fTrackLabel, fOpWaveformLabel, fOpHitLabel, fOpFlashLabel, fGEANTLabel; // --- Input settings imported from the fcl std::vector fLabels, fBackgroundLabels; @@ -88,12 +88,12 @@ namespace solar float fMaxSignalK; float fClusterMatchTime, fAdjClusterRad, fMinClusterCharge, fClusterMatchCharge, fAdjOpFlashX, fAdjOpFlashY, fAdjOpFlashZ, fAdjOpFlashMaxPERatioCut, fAdjOpFlashMinPECut, fAdjOpFlashMinPEAttenuation, fClusterMatchNHit, fClusterAlgoTime; float fOpFlashTimeOffset, fOpFlashAlgoMinTime, fOpFlashAlgoMaxTime, fOpFlashAlgoRad, fOpFlashAlgoPE, fOpFlashAlgoTriggerPE, fOpFlashAlgoHotVertexThld, fXACathodeX, fXAMembraneY, fXAStartCapZ, fXAFinalCapZ; - bool fOpFlashAlgoHitDuplicates; + bool fOpFlashAlgoHitDuplicates, fOpFlashAlgoWeightedTime; bool fClusterPreselectionSignal, fClusterPreselectionPrimary, fClusterPreselectionTrack, fClusterPreselectionFlashMatch; bool fGenerateSolarCluster, fGenerateAdjCluster, fGenerateAdjOpFlash, fFlashMatchByResidual; bool fSaveSignalDaughters, fSaveSignalEDep, fSaveSignalOpHits, fSaveOpFlashInfo, fSaveAdjOpFlashInfo, fSaveTrackInfo; bool fAdjOpFlashMembraneProjection, fAdjOpFlashEndCapProjection; // If true, the TPC reco is projected to the membrane plane. If false, apply a 3D constraint dT, Y, Z. - bool fOpFlashTime2us; // If true, the OpFlash time is in ticks, and we convert it to microseconds. + bool fOpHitTime2us, fOpFlashTime2us; // If true, the OpHit / OpFlash time is in ticks, and we convert it to microseconds. std::vector SelectedEvents; // List of events that pass the selection. 0 = not selected, 1 = selected. // --- Our TTrees, and its associated variables. TTree *fConfigTree; @@ -119,12 +119,15 @@ namespace solar // --- OpFlash Variables std::vector OpFlashID, OpFlashNHits, OpFlashPlane; - std::vector OpFlashPur, OpFlashPE, OpFlashMaxPE, OpFlashX, OpFlashY, OpFlashZ, OpFlashTime, OpFlashDeltaT, OpFlashSTD, OpFlashFast; + std::vector OpHitAmplitude, OpFlashPur, OpFlashPE, OpFlashMaxPE, OpFlashX, OpFlashY, OpFlashZ, OpFlashTime, OpFlashDeltaT, OpFlashSTD, OpFlashFast, OpFlashWaveformTime; + std::vector OpFlashWaveformValid; + std::vector> OpFlashWaveform; // --- MatchedFlash Variables int MFlashNHits, MFlashPlane; - float MFlashR, MFlashPE, MFlashMaxPE, MFlashPur, MFlashFast, MFlashTime, MFlashSTD, MFlashRecoX, MFlashRecoY, MFlashRecoZ, MFlashResidual; - bool MFlashCorrect; + float MOpHitAmplitude, MFlashR, MFlashPE, MFlashMaxPE, MFlashPur, MFlashFast, MFlashTime, MFlashSTD, MFlashRecoX, MFlashRecoY, MFlashRecoZ, MFlashResidual, MFlashWaveformTime; + std::vector MFlashWaveform; + bool MFlashCorrect, MFlashWaveformValid; // --- Maps to hold the geo::TPCID object for each TPCid std::map TPCIDMap; // Key is the TPC index, value is the TPCID object @@ -171,13 +174,15 @@ namespace solar fSolarClusterLabel = p.get("SolarClusterLabel", "solarcluster"); fClusterChargeVariable = p.get("ClusterChargeVariable", "Integral"); fTrackLabel = p.get("TrackLabel", "pmtrack"); + fOpWaveformLabel = p.get("OpWaveformLabel", "opdec"); fOpHitLabel = p.get("OpHitLabel", "ophitspe"); fOpHitTimeVariable = p.get("OpHitTimeVariable", "PeakTime"); fOpFlashLabel = p.get("OpFlashLabel", "solarflash"); - fOpFlashTime2us = p.get("OpFlashTime2us", false); // If true, the OpFlash time is in ticks, and we convert it to microseconds. - fOpFlashTimeOffset = p.get("OpFlashTimeOffset", 0.0); // Time offset to be applied to the OpFlash time (in us if fOpFlashTime2us is true, in ticks otherwise) - fMaxSignalK = p.get("MaxSignalK", 30.0); // Maximum kinetic energy in [MeV] of the particle considered as signal - fClusterAlgoTime = p.get("ClusterAlgoTime", 12.5); // Time window (in us) to look for hits to be clustered together + fOpHitTime2us = p.get("OpHitTime2us", false); // If true, the OpHit time is in ticks, and we convert it to microseconds. + fOpFlashTime2us = p.get("OpFlashTime2us", false); // If true, the OpFlash time is in ticks, and we convert it to microseconds. + fOpFlashTimeOffset = p.get("OpFlashTimeOffset", 0.0); // Time offset to be applied to the OpFlash time (in us if fOpFlashTime2us is true, in ticks otherwise) + fMaxSignalK = p.get("MaxSignalK", 30.0); // Maximum kinetic energy in [MeV] of the particle considered as signal + fClusterAlgoTime = p.get("ClusterAlgoTime", 12.5); // Time window (in us) to look for hits to be clustered together fClusterAlgoAdjChannel = p.get("ClusterAlgoAdjChannel"); fGenerateSolarCluster = p.get("GenerateSolarCluster",true); fClusterMatchNHit = p.get("ClusterMatchNHit", 2.0); @@ -198,8 +203,9 @@ namespace solar fXAMembraneY = p.get("XAMembraneY"); fXAStartCapZ = p.get("XAStartCapZ"); fXAFinalCapZ = p.get("XAFinalCapZ"); - fOpFlashAlgoMinTime = p.get("OpFlashAlgoMinTime", 0.010); // 10 ns [0.6 tick] - fOpFlashAlgoMaxTime = p.get("OpFlashAlgoMaxTime", 0.016); // 16 ns [1 tick] + fOpFlashAlgoMinTime = p.get("OpFlashAlgoMinTime", 0.32); // [20 PDS tick] + fOpFlashAlgoMaxTime = p.get("OpFlashAlgoMaxTime", 0.96); // [60 PDS tick] + fOpFlashAlgoWeightedTime = p.get("OpFlashAlgoWeightedTime"); fOpFlashAlgoRad = p.get("OpFlashAlgoRad"); fOpFlashAlgoPE = p.get("OpFlashAlgoPE"); fOpFlashAlgoTriggerPE = p.get("OpFlashAlgoTriggerPE"); @@ -207,9 +213,9 @@ namespace solar fOpFlashAlgoHitDuplicates = p.get("OpFlashAlgoHitDuplicates"); fAdjOpFlashMembraneProjection = p.get("AdjOpFlashMembraneProjection"); fAdjOpFlashEndCapProjection = p.get("AdjOpFlashEndCapProjection"); - fAdjOpFlashX = p.get("AdjOpFlashX", 100.0); - fAdjOpFlashY = p.get("AdjOpFlashY", 100.0); - fAdjOpFlashZ = p.get("AdjOpFlashZ", 100.0); + fAdjOpFlashX = p.get("AdjOpFlashX", 140.0); + fAdjOpFlashY = p.get("AdjOpFlashY", 140.0); + fAdjOpFlashZ = p.get("AdjOpFlashZ", 140.0); fAdjOpFlashMaxPERatioCut = p.get("AdjOpFlashMaxPERatioCut"); fAdjOpFlashMinPECut = p.get("AdjOpFlashMinPECut"); fAdjOpFlashMinPEAttenuate = p.get("AdjOpFlashMinPEAttenuate"); @@ -250,6 +256,7 @@ namespace solar fConfigTree->Branch("TrackLabel", &fTrackLabel); fConfigTree->Branch("OpHitLabel", &fOpHitLabel); fConfigTree->Branch("OpFlashLabel", &fOpFlashLabel); + fConfigTree->Branch("OpHitTime2us", &fOpHitTime2us); fConfigTree->Branch("OpFlashTime2us", &fOpFlashTime2us); fConfigTree->Branch("OpFlashTimeOffset", &fOpFlashTimeOffset); fConfigTree->Branch("ClusterAlgoTime", &fClusterAlgoTime); @@ -274,6 +281,7 @@ namespace solar fConfigTree->Branch("XAFinalCapZ", &fXAFinalCapZ); fConfigTree->Branch("OpFlashAlgoMinTime", &fOpFlashAlgoMinTime); fConfigTree->Branch("OpFlashAlgoMaxTime", &fOpFlashAlgoMaxTime); + fConfigTree->Branch("OpFlashAlgoWeightedTime", &fOpFlashAlgoWeightedTime); fConfigTree->Branch("OpFlashAlgoRad", &fOpFlashAlgoRad); fConfigTree->Branch("OpFlashAlgoPE", &fOpFlashAlgoPE); fConfigTree->Branch("OpFlashAlgoTriggerPE", &fOpFlashAlgoTriggerPE); @@ -311,7 +319,7 @@ namespace solar fMCTruthTree->Branch("SignalParticleY", &SignalParticleY, "SignalParticleY/F"); // True signal Y [cm] fMCTruthTree->Branch("SignalParticleZ", &SignalParticleZ, "SignalParticleZ/F"); // True signal Z [cm] fMCTruthTree->Branch("SignalParticlePDG", &SignalParticlePDG, "SignalParticlePDG/I"); // True signal PDG - fMCTruthTree->Branch("SignalParticleTime", &SignalParticleTime, "SignalParticleTime/F"); // True signal time [tick] + fMCTruthTree->Branch("SignalParticleTime", &SignalParticleTime, "SignalParticleTime/F"); // True signal time [us] fMCTruthTree->Branch("OpHitNum", &OpHitNum, "OpHitNum/I"); // Number of OpHits fMCTruthTree->Branch("OpFlashNum", &OpFlashNum, "OpFlashNum/I"); // Number of OpFlashes fMCTruthTree->Branch("HitNum", &HitNum, "HitNum/I"); // Number of hits in each TPC plane @@ -324,7 +332,7 @@ namespace solar fMCTruthTree->Branch("TSignalE", &SignalEList); // Energy of Signal particles [MeV] fMCTruthTree->Branch("TSignalP", &SignalPList); // Energy of Signal momentum [MeV] fMCTruthTree->Branch("TSignalK", &SignalKList); // Kinetik Energy of Signal particles [MeV] - fMCTruthTree->Branch("TSignalT", &SignalTimeList); // Time of Signal particles [ticks] + fMCTruthTree->Branch("TSignalT", &SignalTimeList); // Time of Signal particles [us] fMCTruthTree->Branch("TSignalEndX", &SignalEndXList); // X of Signal particles [cm] fMCTruthTree->Branch("TSignalEndY", &SignalEndYList); // Y of Signal particles [cm] fMCTruthTree->Branch("TSignalEndZ", &SignalEndZList); // Z of Signal particles [cm] @@ -370,6 +378,8 @@ namespace solar fMCTruthTree->Branch("OpFlashNHits", &OpFlashNHits); // OpFlash NHit fMCTruthTree->Branch("OpFlashPlane", &OpFlashPlane); // OpFlash Plane fMCTruthTree->Branch("OpFlashMaxPE", &OpFlashMaxPE); // OpFlash Max PE + fMCTruthTree->Branch("OpFlashWaveform", &OpFlashWaveform); // OpFlash Waveform + fMCTruthTree->Branch("OpFlashWaveformValid", &OpFlashWaveformValid); // OpFlash Waveform Valid } // Repeated Truth info. @@ -384,7 +394,7 @@ namespace solar fSolarNuAnaTree->Branch("SignalParticleY", &SignalParticleY, "SignalParticleY/F"); // True signal Y [cm] fSolarNuAnaTree->Branch("SignalParticleZ", &SignalParticleZ, "SignalParticleZ/F"); // True signal Z [cm] fSolarNuAnaTree->Branch("SignalParticlePDG", &SignalParticlePDG, "SignalParticlePDG/I"); // True signal PDG - fSolarNuAnaTree->Branch("SignalParticleTime", &SignalParticleTime, "SignalParticleTime/F"); // True signal Time [tick] + fSolarNuAnaTree->Branch("SignalParticleTime", &SignalParticleTime, "SignalParticleTime/F"); // True signal Time [us] if (fSaveSignalDaughters) { // Save Signal Daughters. (Only makes sense for marley) @@ -411,7 +421,7 @@ namespace solar fSolarNuAnaTree->Branch("Generator", &MGen, "Generator/I"); // Main cluster generator idx fSolarNuAnaTree->Branch("GenPurity", &MGenPur, "GenPurity/F"); // Main cluster reco generator purity fSolarNuAnaTree->Branch("TPC", &MTPC, "ColTPC/I"); // Main cluster TPC - fSolarNuAnaTree->Branch("Time", &MTime, "ColTime/F"); // Main cluster time [ticks] + fSolarNuAnaTree->Branch("Time", &MTime, "ColTime/F"); // Main cluster time [us] fSolarNuAnaTree->Branch("NHits", &MNHit, "ColNHits/I"); // Main cluster #hits fSolarNuAnaTree->Branch("Charge", &MCharge, "ColCharge/F"); // Main cluster charge [ADC*ticks] fSolarNuAnaTree->Branch("MaxCharge", &MMaxCharge, "ColCharge/F"); // Main cluster's max TPCHit-charge [ADC*ticks] @@ -434,13 +444,13 @@ namespace solar fSolarNuAnaTree->Branch("MainE", &MMainE, "MainE/F"); // Main cluster main energy [MeV] fSolarNuAnaTree->Branch("MainP", &MMainP, "MainP/F"); // Main cluster main momentum [MeV] fSolarNuAnaTree->Branch("MainK", &MMainK, "MainK/F"); // Main cluster main kinetic energy [MeV] - fSolarNuAnaTree->Branch("MainTime", &MMainTime, "MainTime/F"); // Main cluster main Time [ticks] + fSolarNuAnaTree->Branch("MainTime", &MMainTime, "MainTime/F"); // Main cluster main Time [us] fSolarNuAnaTree->Branch("MainPDG", &MMainPDG, "MainPDG/I"); // Main cluster main pdg fSolarNuAnaTree->Branch("MainParentPDG", &MMainParentPDG, "MainParentPDG/I"); // Main cluster main pdg fSolarNuAnaTree->Branch("MainParentE", &MMainParentE, "MainParentE/F"); // Main cluster main parent energy [MeV] fSolarNuAnaTree->Branch("MainParentP", &MMainParentP, "MainParentP/F"); // Main cluster main parent momentum [MeV] fSolarNuAnaTree->Branch("MainParentK", &MMainParentK, "MainParentK/F"); // Main cluster main parent kinetic energy [MeV] - fSolarNuAnaTree->Branch("MainParentTime", &MMainParentTime, "MainParentTime/F"); // Main cluster main parent Time [ticks] + fSolarNuAnaTree->Branch("MainParentTime", &MMainParentTime, "MainParentTime/F"); // Main cluster main parent Time [us] fSolarNuAnaTree->Branch("MainVertex", &MMainVertex); // Main cluster main particle vertex [cm] fSolarNuAnaTree->Branch("EndVertex", &MEndVertex); // Main cluster end particle vertex [cm] fSolarNuAnaTree->Branch("MainParentVertex", &MMainParentVertex); // Main cluster parent particle vertex [cm] @@ -461,7 +471,7 @@ namespace solar fSolarNuAnaTree->Branch("AdjClNHits", &MAdjClNHits); // Adj. clusters' #hits fSolarNuAnaTree->Branch("AdjClInd0NHits", &MAdjClInd0NHits); // Adj. clusters' #hits fSolarNuAnaTree->Branch("AdjClInd1NHits", &MAdjClInd1NHits); // Adj. clusters' #hits - fSolarNuAnaTree->Branch("AdjClTime", &MAdjClTime); // Adj. clusters' time [ticks] + fSolarNuAnaTree->Branch("AdjClTime", &MAdjClTime); // Adj. clusters' time [us] fSolarNuAnaTree->Branch("AdjClCharge", &MAdjClCharge); // Adj. clusters' charge [ADC*ticks] fSolarNuAnaTree->Branch("AdjClInd0Charge", &MAdjClInd0Charge); // Adj. clusters' charge [ADC*ticks] fSolarNuAnaTree->Branch("AdjClInd1Charge", &MAdjClInd1Charge); // Adj. clusters' charge [ADC*ticks] @@ -487,36 +497,40 @@ namespace solar // Adj. Flash info. if (fSaveAdjOpFlashInfo) { - fSolarNuAnaTree->Branch("AdjOpFlashR", &MAdjFlashR); // Adj. flash' reco distance [cm] - fSolarNuAnaTree->Branch("AdjOpFlashPE", &MAdjFlashPE); // Adj. flash' tot #PE [ADC*ticks] - fSolarNuAnaTree->Branch("AdjOpFlashPur", &MAdjFlashPur); // Adj. flash' purity - fSolarNuAnaTree->Branch("AdjOpFlashSTD", &MAdjFlashSTD); // Adj. flash' STD - fSolarNuAnaTree->Branch("AdjOpFlashFast", &MAdjFlashFast); // Adj. flash' Fast Component - fSolarNuAnaTree->Branch("AdjOpFlashTime", &MAdjFlashTime); // Adj. flash' time [ticks] - fSolarNuAnaTree->Branch("AdjOpFlashNHits", &MAdjFlashNHits); // Adj. flash' #hits - fSolarNuAnaTree->Branch("AdjOpFlashPlane", &MAdjFlashPlane); // Adj. flash' Plane - fSolarNuAnaTree->Branch("AdjOpFlashMaxPE", &MAdjFlashMaxPE); // Adj. flash' max #PE [ADC*ticks] - fSolarNuAnaTree->Branch("AdjOpFlashRecoX", &MAdjFlashRecoX); // Adj. flash' reco X [cm] - fSolarNuAnaTree->Branch("AdjOpFlashRecoY", &MAdjFlashRecoY); // Adj. flash' reco Y [cm] - fSolarNuAnaTree->Branch("AdjOpFlashRecoZ", &MAdjFlashRecoZ); // Adj. flash' reco Z [cm] - fSolarNuAnaTree->Branch("AdjOpFlashResidual", &MAdjFlashResidual); // Adj. flash' residual wrt. cluster + fSolarNuAnaTree->Branch("AdjOpFlashR", &MAdjFlashR); // Adj. flash reco distance [cm] + fSolarNuAnaTree->Branch("AdjOpFlashPE", &MAdjFlashPE); // Adj. flash tot #PE [PE] + fSolarNuAnaTree->Branch("AdjOpFlashPur", &MAdjFlashPur); // Adj. flash purity + fSolarNuAnaTree->Branch("AdjOpFlashSTD", &MAdjFlashSTD); // Adj. flash STD + fSolarNuAnaTree->Branch("AdjOpFlashFast", &MAdjFlashFast); // Adj. flash Fast Component + fSolarNuAnaTree->Branch("AdjOpFlashTime", &MAdjFlashTime); // Adj. flash time [us] + fSolarNuAnaTree->Branch("AdjOpFlashNHits", &MAdjFlashNHits); // Adj. flash #hits + fSolarNuAnaTree->Branch("AdjOpFlashPlane", &MAdjFlashPlane); // Adj. flash Plane + fSolarNuAnaTree->Branch("AdjOpFlashMaxPE", &MAdjFlashMaxPE); // Adj. flash max #PE [PE] + fSolarNuAnaTree->Branch("AdjOpFlashRecoX", &MAdjFlashRecoX); // Adj. flash reco X [cm] + fSolarNuAnaTree->Branch("AdjOpFlashRecoY", &MAdjFlashRecoY); // Adj. flash reco Y [cm] + fSolarNuAnaTree->Branch("AdjOpFlashRecoZ", &MAdjFlashRecoZ); // Adj. flash reco Z [cm] + fSolarNuAnaTree->Branch("AdjOpFlashResidual", &MAdjFlashResidual); // Adj. flash residual wrt. cluster } // Matched Flash info. - fSolarNuAnaTree->Branch("MatchedOpFlashR", &MFlashR, "MatchedOpFlashR/F"); // Matched flash' reco distance [cm] - fSolarNuAnaTree->Branch("MatchedOpFlashPE", &MFlashPE, "MatchedOpFlashPE/F"); // Matched flash' tot #PE [ADC*ticks] - fSolarNuAnaTree->Branch("MatchedOpFlashPur", &MFlashPur, "MatchedOpFlashPur/F"); // Matched flash' purity - fSolarNuAnaTree->Branch("MatchedOpFlashSTD", &MFlashSTD, "MatchedOpFlashSTD/F"); // Matched flash' STD - fSolarNuAnaTree->Branch("MatchedOpFlashFast", &MFlashFast, "MatchedOpFlashFast/F"); // Matched flash' Fast Component - fSolarNuAnaTree->Branch("MatchedOpFlashTime", &MFlashTime, "MatchedOpFlashTime/F"); // Matched flash' time [ticks] - fSolarNuAnaTree->Branch("MatchedOpFlashNHits", &MFlashNHits, "MatchedOpFlashNHits/I"); // Matched flash' #hits - fSolarNuAnaTree->Branch("MatchedOpFlashPlane", &MFlashPlane, "MatchedOpFlashPlane/I"); // Matched flash' Plane - fSolarNuAnaTree->Branch("MatchedOpFlashMaxPE", &MFlashMaxPE, "MatchedOpFlashMaxPE/F"); // Matched flash' max #PE [ADC*ticks] - fSolarNuAnaTree->Branch("MatchedOpFlashRecoX", &MFlashRecoX, "MatchedOpFlashRecoX/F"); // Matched flash' reco X [cm] - fSolarNuAnaTree->Branch("MatchedOpFlashRecoY", &MFlashRecoY, "MatchedOpFlashRecoY/F"); // Matched flash' reco Y [cm] - fSolarNuAnaTree->Branch("MatchedOpFlashRecoZ", &MFlashRecoZ, "MatchedOpFlashRecoZ/F"); // Matched flash' reco Z [cm] - fSolarNuAnaTree->Branch("MatchedOpFlashResidual", &MFlashResidual, "MatchedOpFlashResidual/F"); // Matched flash' residual wrt. cluster - fSolarNuAnaTree->Branch("MatchedOpFlashCorrectly", &MFlashCorrect); // Matched flash' correctnes (bool) + fSolarNuAnaTree->Branch("MatchedOpHitAmplitude", &MOpHitAmplitude, "MatchedOpHitAmplitude/F"); // Matched OpHit amplitude [ADC] for raw [PE] for deconvolved wvfs + fSolarNuAnaTree->Branch("MatchedOpFlashR", &MFlashR, "MatchedOpFlashR/F"); // Matched flash reco distance [cm] + fSolarNuAnaTree->Branch("MatchedOpFlashPE", &MFlashPE, "MatchedOpFlashPE/F"); // Matched flash tot #PE [PE] + fSolarNuAnaTree->Branch("MatchedOpFlashPur", &MFlashPur, "MatchedOpFlashPur/F"); // Matched flash purity + fSolarNuAnaTree->Branch("MatchedOpFlashSTD", &MFlashSTD, "MatchedOpFlashSTD/F"); // Matched flash STD + fSolarNuAnaTree->Branch("MatchedOpFlashFast", &MFlashFast, "MatchedOpFlashFast/F"); // Matched flash Fast Component + fSolarNuAnaTree->Branch("MatchedOpFlashTime", &MFlashTime, "MatchedOpFlashTime/F"); // Matched flash time [us] + fSolarNuAnaTree->Branch("MatchedOpFlashNHits", &MFlashNHits, "MatchedOpFlashNHits/I"); // Matched flash #hits + fSolarNuAnaTree->Branch("MatchedOpFlashPlane", &MFlashPlane, "MatchedOpFlashPlane/I"); // Matched flash Plane + fSolarNuAnaTree->Branch("MatchedOpFlashMaxPE", &MFlashMaxPE, "MatchedOpFlashMaxPE/F"); // Matched flash max #PE [PE] + fSolarNuAnaTree->Branch("MatchedOpFlashRecoX", &MFlashRecoX, "MatchedOpFlashRecoX/F"); // Matched flash reco X [cm] + fSolarNuAnaTree->Branch("MatchedOpFlashRecoY", &MFlashRecoY, "MatchedOpFlashRecoY/F"); // Matched flash reco Y [cm] + fSolarNuAnaTree->Branch("MatchedOpFlashRecoZ", &MFlashRecoZ, "MatchedOpFlashRecoZ/F"); // Matched flash reco Z [cm] + fSolarNuAnaTree->Branch("MatchedOpFlashResidual", &MFlashResidual, "MatchedOpFlashResidual/F"); // Matched flash residual wrt. cluster + fSolarNuAnaTree->Branch("MatchedOpFlashWaveform", &MFlashWaveform); // Matched flash waveform + fSolarNuAnaTree->Branch("MatchedOpFlashWaveformTime", &MFlashWaveformTime); // Matched flash waveform time [us] + fSolarNuAnaTree->Branch("MatchedOpFlashWaveformValid", &MFlashWaveformValid); // Matched flash waveform valid + fSolarNuAnaTree->Branch("MatchedOpFlashCorrectly", &MFlashCorrect); // Matched flash correctnes (bool) fConfigTree->AddFriend(fSolarNuAnaTree); fMCTruthTree->AddFriend(fSolarNuAnaTree); @@ -893,6 +907,8 @@ namespace solar //------------------------------------------------------------------- Optical Flash Analysis --------------------------------------------------------------------// //---------------------------------------------------------------------------------------------------------------------------------------------------------------// // Find OpHits and OpFlashes associated with the event + float MinOpHitTime = 1e6, MinOpFlashTime = 1e6; + float MaxOpHitTime = -1e6, MaxOpFlashTime = -1e6; std::string sOpFlashTruth = ""; std::vector> OpHitList; art::Handle> OpHitHandle; @@ -909,55 +925,78 @@ namespace solar adjophits->CalcAdjOpHits(OpHitList, OpHitVec, OpHitIdx, evt); adjophits->MakeFlashVector(FlashVec, OpHitVec, evt); OpFlashNum = int(FlashVec.size()); - for (int i = 0; i < int(FlashVec.size()); i++) { AdjOpHitsUtils::FlashInfo TheFlash = FlashVec[i]; double ThisOpFlashPur = 0; + float ThisOpFlashTime = -1e6; + + if (fOpFlashAlgoWeightedTime) + ThisOpFlashTime = TheFlash.TimeWeighted - fOpFlashTimeOffset; + else + ThisOpFlashTime = TheFlash.MainOpHitTime - fOpFlashTimeOffset; + + OpHitAmplitude.push_back(TheFlash.MainOpHitAmplitude); OpFlashPlane.push_back(TheFlash.Plane); OpFlashNHits.push_back(TheFlash.NHit); - OpFlashTime.push_back(TheFlash.Time - fOpFlashTimeOffset); // Convert to microseconds happens in AdjOpHits + OpFlashTime.push_back(ThisOpFlashTime); OpFlashDeltaT.push_back(TheFlash.TimeWidth); // Convert to microseconds OpFlashPE.push_back(TheFlash.PE); - OpFlashMaxPE.push_back(TheFlash.MaxPE); + OpFlashMaxPE.push_back(TheFlash.MainOpHitPE); OpFlashFast.push_back(TheFlash.FastToTotal); OpFlashID.push_back(i); OpFlashX.push_back(TheFlash.X); OpFlashY.push_back(TheFlash.Y); OpFlashZ.push_back(TheFlash.Z); OpFlashSTD.push_back(TheFlash.STD); + OpFlashWaveform.push_back(TheFlash.MainOpWaveform); + OpFlashWaveformTime.push_back(TheFlash.MainOpWaveformTime); + OpFlashWaveformValid.push_back(TheFlash.MainOpWaveformValid); + for (int j = 0; j < int(OpHitVec[i].size()); j++) { recob::OpHit OpHit = *OpHitVec[i][j]; + const std::vector ThisOpHitTrackIds = pbt->OpHitToTrackIds(OpHit); - float ThisOphitPurity = 0; + float ThisOpHitTime = -1e6; + float ThisOpHitPurity = 0; for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) { if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) - ThisOphitPurity += 1; + ThisOpHitPurity += 1; } // Check if ThisOpHitTrackIds is empty if (ThisOpHitTrackIds.size() == 0) - ThisOphitPurity = 0; + ThisOpHitPurity = 0; else - ThisOphitPurity /= int(ThisOpHitTrackIds.size()); + ThisOpHitPurity /= int(ThisOpHitTrackIds.size()); - ThisOpFlashPur += ThisOphitPurity * OpHit.PE(); + ThisOpFlashPur += ThisOpHitPurity * OpHit.PE(); auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(OpHit.OpChannel()).GetCenter(); - SOpHitPur.push_back(ThisOphitPurity); + SOpHitPur.push_back(ThisOpHitPurity); SOpHitChannel.push_back(OpHit.OpChannel()); if (fOpHitTimeVariable == "StartTime") - SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + ThisOpHitTime = OpHit.StartTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds else // Default to PeakTime - SOpHitTime.push_back(OpHit.PeakTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds - + ThisOpHitTime = OpHit.PeakTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.PeakTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + if (fOpHitTime2us) + ThisOpHitTime *= clockData.OpticalClock().TickPeriod(); // Convert to microseconds + + ThisOpHitTime -= fOpFlashTimeOffset; + SOpHitTime.push_back(ThisOpHitTime); SOpHitPE.push_back(OpHit.PE()); SOpHitX.push_back(OpHitXYZ.X()); SOpHitY.push_back(OpHitXYZ.Y()); SOpHitZ.push_back(OpHitXYZ.Z()); SOpHitFlashID.push_back(i); SOpHitPlane.push_back(TheFlash.Plane); + if (ThisOpHitTime < MinOpHitTime) + MinOpHitTime = ThisOpHitTime; + if (ThisOpHitTime > MaxOpHitTime) + MaxOpHitTime = ThisOpHitTime; } // Check if OpHitVec[i] is empty if (OpHitVec[i].size() == 0) @@ -966,15 +1005,26 @@ namespace solar ThisOpFlashPur /= TheFlash.PE; OpFlashPur.push_back(ThisOpFlashPur); + if (ThisOpFlashTime < MinOpFlashTime) + MinOpFlashTime = ThisOpFlashTime; + if (ThisOpFlashTime > MaxOpFlashTime) + MaxOpFlashTime = ThisOpFlashTime; + if (ThisOpFlashPur > 0) { - sOpFlashTruth += "OpFlash PE " + ProducerUtils::str(TheFlash.PE) + " with purity " + ProducerUtils::str(ThisOpFlashPur) + " time " + ProducerUtils::str(TheFlash.Time) + " plane " + ProducerUtils::str(TheFlash.Plane) + "\n"; + sOpFlashTruth += "OpFlash PE " + ProducerUtils::str(TheFlash.PE) + " with purity " + ProducerUtils::str(ThisOpFlashPur) + " time " + ProducerUtils::str(ThisOpFlashTime) + " plane " + ProducerUtils::str(TheFlash.Plane) + "\n"; sOpFlashTruth += " - Vertex (" + ProducerUtils::str(TheFlash.X) + ", " + ProducerUtils::str(TheFlash.Y) + ", " + ProducerUtils::str(TheFlash.Z) + ")\n"; - sOpFlashTruth += "\t*** 1st Sanity check: Ratio " + ProducerUtils::str(TheFlash.MaxPE / TheFlash.PE) + " <= " + ProducerUtils::str(fAdjOpFlashMaxPERatioCut) + "\n"; + sOpFlashTruth += "\t*** 1st Sanity check: Ratio " + ProducerUtils::str(TheFlash.MainOpHitPE / TheFlash.PE) + " <= " + ProducerUtils::str(fAdjOpFlashMaxPERatioCut) + "\n"; sOpFlashTruth += "\t*** 2nd Sanity check: #OpHits " + ProducerUtils::str(int(OpHitVec[i].size())) + " >= " + ProducerUtils::str(TheFlash.NHit) + "\n"; } } } else { + float MainWaveformTime = -1e6; + bool MainWaveformValid = false; + std::vector MainWaveform; + std::vector MatchedValidWaveforms; + std::vector MatchedWaveformTimes; + std::vector> MatchedWaveforms; std::vector> OpFlashList; art::Handle> FlashHandle; @@ -988,7 +1038,10 @@ namespace solar for (int i = 0; i < int(OpFlashList.size()); i++) { recob::OpFlash TheFlash = *OpFlashList[i]; + std::vector> MatchedHits = OpAssns.at(i); + adjophits->GetOpHitSignal(MatchedHits, MatchedWaveforms, MatchedWaveformTimes, MatchedValidWaveforms, evt); + int NMatchedHits = MatchedHits.size(); double FlashStdDev = 0.0, TotalFlashPE = 0, MaxOpHitPE = 0; std::vector varXY, varYZ, varXZ; @@ -1000,12 +1053,13 @@ namespace solar art::Ptr OpHitPtr = MatchedHits[j]; mf::LogDebug("SolarNuAna") << "Assigning OpHit to Flash"; const std::vector ThisOpHitTrackIds = pbt->OpHitToTrackIds(OpHit); - float ThisOphitPurity = 0; + float ThisOpHitPurity = 0; + float ThisOpHitTime = -1e6; for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) { if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) { - ThisOphitPurity += 1; + ThisOpHitPurity += 1; } } @@ -1015,10 +1069,20 @@ namespace solar varXY.push_back(sqrt(pow(TheFlash.XCenter() - OpHitXYZ.X(), 2) + pow(TheFlash.YCenter() - OpHitXYZ.Y(), 2)) * OpHit.PE()); varYZ.push_back(sqrt(pow(TheFlash.YCenter() - OpHitXYZ.Y(), 2) + pow(TheFlash.ZCenter() - OpHitXYZ.Z(), 2)) * OpHit.PE()); varXZ.push_back(sqrt(pow(TheFlash.XCenter() - OpHitXYZ.X(), 2) + pow(TheFlash.ZCenter() - OpHitXYZ.Z(), 2)) * OpHit.PE()); - SOpHitPur.push_back(ThisOphitPurity / int(ThisOpHitTrackIds.size())); + SOpHitPur.push_back(ThisOpHitPurity / int(ThisOpHitTrackIds.size())); if (OpHit.PE() > MaxOpHitPE) { MaxOpHitPE = OpHit.PE(); + if (MatchedValidWaveforms[j]) { + MainWaveform = MatchedWaveforms[j]; + MainWaveformTime = MatchedWaveformTimes[j]; + MainWaveformValid = MatchedValidWaveforms[j]; + } + else { + MainWaveform = {}; + MainWaveformTime = -1e6; + MainWaveformValid = false; + } } SOpHitFlashID.push_back(i); @@ -1028,12 +1092,22 @@ namespace solar SOpHitZ.push_back(OpHitXYZ.Z()); if (fOpHitTimeVariable == "StartTime") - SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + ThisOpHitTime = OpHit.StartTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds else // Default to PeakTime - SOpHitTime.push_back(OpHit.PeakTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds - + ThisOpHitTime = OpHit.PeakTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.PeakTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + if (fOpHitTime2us) + ThisOpHitTime *= clockData.OpticalClock().TickPeriod(); // Convert to microseconds + + ThisOpHitTime -= fOpFlashTimeOffset; SOpHitChannel.push_back(OpHit.OpChannel()); SOpHitPlane.push_back(adjophits->GetOpHitPlane(OpHitPtr, 0.1)); // Get plane assignment for the OpHit + + if (ThisOpHitTime < MinOpHitTime) + MinOpHitTime = ThisOpHitTime; + if (ThisOpHitTime > MaxOpHitTime) + MaxOpHitTime = ThisOpHitTime; } // End of OpHit loop OpHitVec.push_back(MatchedHits); @@ -1054,16 +1128,24 @@ namespace solar OpFlashPE.push_back(TheFlash.TotalPE()); OpFlashNHits.push_back(MatchedHits.size()); OpFlashFast.push_back(TheFlash.FastToTotal()); + OpFlashWaveform.push_back(MainWaveform); + OpFlashWaveformTime.push_back(MainWaveformTime); + OpFlashWaveformValid.push_back(MainWaveformValid); if (fOpFlashTime2us) { OpFlashTime.push_back(TheFlash.Time() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Expected flash to provide time in ticks, convert to microseconds OpFlashDeltaT.push_back(TheFlash.TimeWidth() * clockData.OpticalClock().TickPeriod()); // Expected flash to provide time width in ticks, convert to microseconds } else { - OpFlashTime.push_back(TheFlash.Time()); // Expected flash to provide time in microseconds + OpFlashTime.push_back(TheFlash.Time() - fOpFlashTimeOffset); // Expected flash to provide time in microseconds OpFlashDeltaT.push_back(TheFlash.TimeWidth()); // Expected flash to provide time width in microseconds } + if (TheFlash.Time() < MinOpFlashTime) + MinOpFlashTime = TheFlash.Time(); + if (TheFlash.Time() > MaxOpFlashTime) + MaxOpFlashTime = TheFlash.Time(); + if (ThisOpFlashPur > 0) { mf::LogDebug("SolarNuAna") << "OpFlash PE " << TheFlash.TotalPE() << " with purity " << ThisOpFlashPur << " time " << TheFlash.Time(); sOpFlashTruth += "OpFlash PE " + ProducerUtils::str(TheFlash.TotalPE()) + " with purity " + ProducerUtils::str(ThisOpFlashPur) + " time " + ProducerUtils::str(TheFlash.Time()) + " plane " + ProducerUtils::str(int(TheFlash.Frame())) + "\n"; @@ -1073,7 +1155,8 @@ namespace solar } } } - sOpFlashTruth = sOpFlashTruth + "\n# of OpFlashes (" + fOpFlashLabel + ") in full geometry: " + ProducerUtils::str(OpFlashNum) + "\n"; + sOpFlashTruth = sOpFlashTruth + "\n# of OpHits (" + fOpHitLabel + "): " + ProducerUtils::str(OpHitNum) + " for times between " + ProducerUtils::str(MinOpHitTime) + " us and " + ProducerUtils::str(MaxOpHitTime) + " us.\n"; + sOpFlashTruth = sOpFlashTruth + "# of OpFlashes (" + fOpFlashLabel + "): " + ProducerUtils::str(OpFlashNum) + " for times between " + ProducerUtils::str(MinOpFlashTime) + " us and " + ProducerUtils::str(MaxOpFlashTime) + " us.\n"; producer->PrintInColor(sOpFlashTruth, ProducerUtils::GetColor("blue")); //---------------------------------------------------------------------------------------------------------------------------------------------------------------// @@ -1159,6 +1242,9 @@ namespace solar AllPlaneClusters = {Clusters0, Clusters1, Clusters2}; } + HitNum = {int(Ind0Hits.size()), int(Ind1Hits.size()), int(ColHits.size())}; + ClusterNum = {int(Clusters0.size()), int(Clusters1.size()), int(Clusters2.size())}; + std::vector>> ClVecGenPur = {{}, {}, {}}; std::vector> ClMainID = {{}, {}, {}}, ClTPC = {{}, {}, {}}, ClNHits = {{}, {}, {}}, ClGen = {{}, {}, {}}; std::vector> ClCharge = {{}, {}, {}}, ClMaxCharge = {{}, {}, {}}, ClT = {{}, {}, {}}, ClX = {{}, {}, {}}, ClY = {{}, {}, {}}, ClZ = {{}, {}, {}}; @@ -1345,7 +1431,6 @@ namespace solar SolarClusterInfo = SolarClusterInfo + "\nFound " + ProducerUtils::str(int(MatchedClustersIdx[2].size())) + " MatchedClusters (from col. plane loop)!"; for (int ThisClIdx = 0; ThisClIdx < int(MatchedClustersIdx[2].size()); ThisClIdx++) { - // MVecTime[2][ThisClIdx] *= clockData.TPCClock().TickPeriod(); // Convert to microseconds for (int plane = 0; plane < 2; plane++) { if (MVecTime[plane][ThisClIdx] > -1e6) { // There is a matched cluster in this plane @@ -1578,7 +1663,7 @@ namespace solar if (fClusterPreselectionPrimary && !MPrimary) { continue; } if (MPrimary) { - sClusterReco += "*** Matched preselection cluster: " + ProducerUtils::str(i) + "\n"; + sClusterReco += "*** Matched preselection cluster: " + ProducerUtils::str(i) + " from " + ProducerUtils::str(int(MVecNHits[2].size())) + "\n"; sClusterReco += " - MainTrackID " + ProducerUtils::str(MVecMainID[2][i]) + "\n"; if (MVecGen[i] > 0 && int(MVecGen[i]) < (int(fLabels.size()) + 1)) { sClusterReco += " - Gen " + ProducerUtils::str(int(MVecGen[i])) + " -> " + fLabels[MVecGen[i] - 1]; @@ -1699,11 +1784,12 @@ namespace solar if ( OpFlashNHits[j] < fAdjOpFlashMinNHitCut || OpFlashMaxPE[j] / OpFlashPE[j] > fAdjOpFlashMaxPERatioCut ) { continue; } - if ( lowe->SelectPDSFlashPE(TPCIDdriftTime[MVecTPC[2][i]], MVecCharge[2][i] - OpFlashTime[j], MVecCharge[2][i], OpFlashPE[j]) == false ) { + if ( lowe->SelectPDSFlashPE(TPCIDdriftTime[MVecTPC[2][i]], MVecTime[2][i] - OpFlashTime[j], MVecCharge[2][i], OpFlashPE[j]) == false ) { continue; // Skip flashes that don't pass the PE selection } // If the residual is smaller than the minimum residual, update the minimum residual and the matched flash. if ((fFlashMatchByResidual && OpFlashResidual < MatchedOpFlashResidual) || (!fFlashMatchByResidual && OpFlashPE[j] > MatchedOpFlashPE)) { + MOpHitAmplitude = OpHitAmplitude[j]; MFlashR = OpFlashR; MFlashPE = OpFlashPE[j]; MFlashFast = OpFlashFast[j]; @@ -1716,17 +1802,22 @@ namespace solar MFlashRecoX = OpFlashX[j]; MFlashRecoY = OpFlashY[j]; MFlashRecoZ = OpFlashZ[j]; + MFlashWaveform = OpFlashWaveform[j]; + MFlashWaveformValid = OpFlashWaveformValid[j]; + MFlashWaveformTime = OpFlashWaveformTime[j]; MFlashResidual = OpFlashResidual; // Create an output string with the flash information. - sFlashReco = "*** Matched flash: \n - Purity " + ProducerUtils::str(OpFlashPur[j]) + + sFlashReco = "*** Matched flash: " + ProducerUtils::str(j) + " from " + ProducerUtils::str(int(OpFlashPE.size())) + "\n" + + " - Purity " + ProducerUtils::str(100*OpFlashPur[j]) + " %" + " Plane " + ProducerUtils::str(OpFlashPlane[j]) + " #Hits " + ProducerUtils::str(OpFlashNHits[j]) + " PE " + ProducerUtils::str(OpFlashPE[j]) + - " MaxPE " + ProducerUtils::str(OpFlashMaxPE[j]) + "\n" + - " - Time " + ProducerUtils::str(OpFlashTime[j]) + - " Fast " + ProducerUtils::str(OpFlashFast[j]) + + " MainOpHitPE " + ProducerUtils::str(OpFlashMaxPE[j]) + "\n" + + " - Time " + ProducerUtils::str(OpFlashTime[j]) + " (us)" + + " Fast " + ProducerUtils::str(100*OpFlashFast[j]) + " %" + " Residual " + ProducerUtils::str(OpFlashResidual) + "\n" + - " - Reco Time,Y,Z ( " + ProducerUtils::str(MFlashTime) + ", " + ProducerUtils::str(OpFlashY[j]) + ", " + ProducerUtils::str(OpFlashZ[j]) + " )" + "\n"; + " - Reco Time,Y,Z ( " + ProducerUtils::str(MFlashTime) + ", " + ProducerUtils::str(OpFlashY[j]) + ", " + ProducerUtils::str(OpFlashZ[j]) + " )" + "\n" + + " - Found valid waveform: " + ProducerUtils::str(MFlashWaveformValid) + "\n"; MatchedOpFlashX = MAdjFlashX; MatchedOpFlashResidual = OpFlashResidual; MatchedOpFlashPE = MFlashPE; @@ -1892,6 +1983,7 @@ namespace solar MFlashRecoY = -1e6; MFlashRecoZ = -1e6; MFlashResidual = -1e6; + MFlashWaveform = {}; MFlashCorrect = false; SignalParticleE = 0; SignalParticleP = 0; @@ -1901,10 +1993,14 @@ namespace solar SignalParticleZ = 0; SignalParticlePDG = 0; SignalParticleTime = 0; + OpHitAmplitude.clear(); OpFlashPur.clear(); OpFlashID.clear(); OpFlashPE.clear(); OpFlashSTD.clear(); + OpFlashWaveform.clear(); + OpFlashWaveformTime.clear(); + OpFlashWaveformValid.clear(); OpFlashFast.clear(); OpFlashMaxPE.clear(); OpFlashX.clear(); diff --git a/duneana/SolarNuAna/fcl/SolarNuAna.fcl b/duneana/SolarNuAna/fcl/SolarNuAna.fcl index 4b5ac724..5348742f 100644 --- a/duneana/SolarNuAna/fcl/SolarNuAna.fcl +++ b/duneana/SolarNuAna/fcl/SolarNuAna.fcl @@ -7,12 +7,14 @@ BEGIN_PROLOG { module_type: "SolarNuAna" GEANT4Label: "largeant" # The label for the process which ran GEANT4 + OpWaveformLabel: "opdec" # The label for the process which ran the OpDetWaveform HitLabel: "hitfd" # String for the process that made the reco hits ClusterLabel: "planecluster" # The label for the process which ran the "perplane" hit clustering SolarClusterLabel: "solarcluster" # The label for the process which ran the solar clustering TrackLabel: "pmtracktc" # The label for the process which ran the PMTrack OpHitLabel: "ophitspe" # The label for the process which ran the OpHit OpFlashLabel: "solarflash" # The label for the process which ran the OpFlash + OpHitTime2us: false # If true, the OpHit time is in [ticks] and will be converted to [us]. OpFlashTime2us: false # If true, the OpFlash time is in [ticks] and will be converted to [us]. For other producers than "solarflash". OpFlashTimeOffset: 18.1 # Time offset to be applied to the OpFlash time in [us] units. @@ -53,18 +55,19 @@ BEGIN_PROLOG GenerateAdjOpFlash: true # Generate OpFlashes. OpFlashAlgoNHit: 0 # Min number of hits to consider a flash. Change to 3 for bkg run to avoid huge output. - OpFlashAlgoMinTime: 0.008 # Negative time window to look for adj. OpHits in [us] units. - OpFlashAlgoMaxTime: 0.016 # Positive time window to look for adj. OpHits in [us] units. - OpFlashAlgoRad: 300 # Distance to look for adj. OpHits in [cm] units. - OpFlashAlgoPE: 1.5 # PE threshold to look for adj. OpHits. - OpFlashAlgoTriggerPE: 1.5 # PE threshold to trigger an OpFlash. + OpFlashAlgoMinTime: 0.32 # Negative time window to look for adj. OpHits in [us] units. + OpFlashAlgoMaxTime: 0.96 # Positive time window to look for adj. OpHits in [us] units. + OpFlashAlgoWeightedTime: true # If true, use weighted time for OpFlash time calculation, otherwise use flash trigger time. + OpFlashAlgoRad: 600 # Distance to look for adj. OpHits in [cm] units. + OpFlashAlgoPE: 0.0 # PE threshold to look for adj. OpHits. + OpFlashAlgoTriggerPE: 0.0 # PE threshold to trigger an OpFlash. OpFlashAlgoHotVertexThld: 0.3 # Relative threshold to consider a hit as hot for opflash vertex determination [0-1]. OpFlashAlgoHitDuplicates: true # If true, allow hits to be used in multiple flashes. AdjOpFlashX: 140. # X distance to search for adj. OpFlashes reconstructed in [cm] units. AdjOpFlashY: 140. # Y distance to search for adj. OpFlashes reconstructed in [cm] units. AdjOpFlashZ: 140. # Z distance to search for adj. OpFlashes reconstructed in [cm] units. - AdjOpFlashMinNHitCut: 3 # Cut on the minimum number of OpHits in the OpFlash. + AdjOpFlashMinNHitCut: 0 # Cut on the minimum number of OpHits in the OpFlash. AdjOpFlashMembraneProjection: true # If true, the OpFlash matching is projected on the membrane planes for VD. AdjOpFlashEndCapProjection: false # If true, the OpFlash matching is projected on the end cap planes for VD. AdjOpFlashMaxPERatioCut: 1.00 # Cut on the maximum OpHit PE contribution to the total OpFlash PE. @@ -103,6 +106,7 @@ BEGIN_PROLOG solar_nu_ana_vd: @local::solar_nu_ana_hd_centralAPA solar_nu_ana_vd.HitLabel: "gaushit" + solar_nu_ana_vd.OpWaveformLabel: "opdigi10ppm" solar_nu_ana_vd.OpHitLabel: "ophit10ppm" solar_nu_ana_vd.OpFlashLabel: "solarflash" solar_nu_ana_vd.OpFlashTimeOffset: 0 @@ -112,8 +116,8 @@ BEGIN_PROLOG solar_nu_ana_vd.XAMembraneY: 743.302 solar_nu_ana_vd.XAFinalCapZ: 2188.38 solar_nu_ana_vd.XAStartCapZ: -96.5 - solar_nu_ana_vd.OpFlashAlgoMinTime: 0.03 - solar_nu_ana_vd.OpFlashAlgoMaxTime: 0.06 + solar_nu_ana_vd.OpFlashAlgoMinTime: 0.96 + solar_nu_ana_vd.OpFlashAlgoMaxTime: 3.20 solar_nu_ana_vd.OpFlashAlgoRad: 600. solar_nu_ana_vd.AdjOpFlashMinPELightMap: [ ["Amplitude", [1.7578e-7, -4.9589e-4, 1.2519e0]], diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl index 08f73118..fc78f1ff 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl @@ -6,6 +6,9 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 +physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 3 physics.analyzers.solarnuana.AdjOpFlashX: 100. physics.analyzers.solarnuana.AdjOpFlashY: 100. physics.analyzers.solarnuana.AdjOpFlashZ: 100. \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl index 47332320..4c2f6ed7 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl @@ -6,6 +6,9 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 +physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 3 physics.analyzers.solarnuana.AdjOpFlashX: 100. physics.analyzers.solarnuana.AdjOpFlashY: 100. physics.analyzers.solarnuana.AdjOpFlashZ: 100. \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl index 8f5ad24e..7a32c5a4 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl @@ -1,11 +1,7 @@ # Specific configuration for investigating energy depositions of signal events. # Run SolarNuAna on the output of a standard reco workflow for DUNE FD -#include "solar_ana_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl" +#include "solar_ana_marley_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl" + +physics.analyzers.solarnuana.SaveOpFlashInfo: false -physics.analyzers.solarnuana.SaveSignalDaughters: true -physics.analyzers.solarnuana.GenerateAdjOpFlash: true -physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true -physics.analyzers.solarnuana.AdjOpFlashX: 100. -physics.analyzers.solarnuana.AdjOpFlashY: 100. -physics.analyzers.solarnuana.AdjOpFlashZ: 100. \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl index e0803a00..518ca233 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_marley_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl @@ -1,11 +1,6 @@ # Specific configuration for investigating energy depositions of signal events. # Run SolarNuAna on the output of a standard reco workflow for DUNE FD -#include "solar_ana_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl" +#include "solar_ana_marley_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl" -physics.analyzers.solarnuana.SaveSignalDaughters: true -physics.analyzers.solarnuana.GenerateAdjOpFlash: true -physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true -physics.analyzers.solarnuana.AdjOpFlashX: 100. -physics.analyzers.solarnuana.AdjOpFlashY: 100. -physics.analyzers.solarnuana.AdjOpFlashZ: 100. \ No newline at end of file +physics.analyzers.solarnuana.SaveOpFlashInfo: false \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_dunevd10kt_1x8x14_3view_30deg.fcl index 01ae7a9c..785851ab 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_dunevd10kt_1x8x14_3view_30deg.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_dunevd10kt_1x8x14_3view_30deg.fcl @@ -5,4 +5,6 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true -physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true \ No newline at end of file +physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl index 8e5dd91e..e3ae4112 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl @@ -6,8 +6,10 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: false -physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 6 +physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 3 physics.analyzers.solarnuana.AdjOpFlashX: 100. physics.analyzers.solarnuana.AdjOpFlashY: 100. physics.analyzers.solarnuana.AdjOpFlashZ: 100. \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl index 36c3c54e..2feee481 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -6,6 +6,8 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: false physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 3 physics.analyzers.solarnuana.AdjOpFlashX: 100. diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl new file mode 100644 index 00000000..b0479de0 --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -0,0 +1,6 @@ +# Specific configuration for investigating energy depositions of signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD-VD + +#include "solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" + +physics.analyzers.solarnuana.SaveSignalDaughters: true \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..7df1b164 --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,6 @@ +# Specific configuration for investigating energy depositions of signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD-VD + +#include "solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.solarnuana.SaveSignalDaughters: true \ No newline at end of file