diff --git a/duneana/MCParticleAna/fcl/MCParticleAna.fcl b/duneana/MCParticleAna/fcl/MCParticleAna.fcl index ef1c9b1a..dd1fbb8b 100644 --- a/duneana/MCParticleAna/fcl/MCParticleAna.fcl +++ b/duneana/MCParticleAna/fcl/MCParticleAna.fcl @@ -26,6 +26,9 @@ BEGIN_PROLOG solar_nu_ana_vd_1x8x14_optimistic: @local::solar_nu_ana_vd_1x8x14 solar_nu_ana_vd_1x8x14_optimistic.ParticleLabelVector: @local::generator_dunevd10kt_1x8x14_3view_30deg_optimistic + solar_nu_ana_vd_1x8x14_shielded: @local::solar_nu_ana_vd_1x8x14 + solar_nu_ana_vd_1x8x14_shielded.ParticleLabelVector: @local::generator_dunevd10kt_1x8x14_3view_30deg_shielded + solar_nu_ana_vd_1x8x6: @local::solar_nu_ana_vd_1x8x14 solar_nu_ana_vd_1x8x6.ParticleLabelVector: @local::generator_dunevd10kt_1x8x6_3view_30deg diff --git a/duneana/MCParticleAna/fcl/mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/MCParticleAna/fcl/mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..082ae9f4 --- /dev/null +++ b/duneana/MCParticleAna/fcl/mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,39 @@ +# mcparticle_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl +# +# Run MCParticleAna on the output of a standard gen file for DUNE FD + +#include "MCParticleAna.fcl" +#include "services_dune.fcl" +#include "tools_dune.fcl" + +process_name: MCParticleAna + +services: +{ + @table::dunefd_services + TFileService: { fileName: "mcparticle_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded_hist.root" } + TimeTracker: {} + MemoryTracker: {} # default is one + RandomNumberGenerator: {} # ART native random number generator + FileCatalogMetadata: @local::art_file_catalog_mc + @table::dunefdvd_1x8x14_3view_30deg_services +} + +physics: +{ + analyzers: + { + mcparticleana: @local::solar_nu_ana_vd_1x8x14_shielded + } + ana: [ mcparticleana ] + end_paths: [ ana ] +} + +source: +{ + module_type: RootInput + maxEvents: -1 # Number of events to create +} + +services.message.destinations.LogStandardOut.threshold: "INFO" +services.message.destinations.LogStandardOut.type: "cout" diff --git a/duneana/MCParticleAna/fcl/mcparticle_ana_alpha_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/MCParticleAna/fcl/mcparticle_ana_alpha_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..ce949800 --- /dev/null +++ b/duneana/MCParticleAna/fcl/mcparticle_ana_alpha_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,4 @@ +#include "mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.mcparticleana.MCParticlePDG: 1000020040 # For alpha particles +physics.analyzers.mcparticleana.MCParticleMinKE: 4. # MeV \ No newline at end of file diff --git a/duneana/MCParticleAna/fcl/mcparticle_ana_electron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/MCParticleAna/fcl/mcparticle_ana_electron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..87682c8f --- /dev/null +++ b/duneana/MCParticleAna/fcl/mcparticle_ana_electron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,4 @@ +#include "mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.mcparticleana.MCParticlePDG: 11 # For electron particles +physics.analyzers.mcparticleana.MCParticleMinKE: 4. # MeV \ No newline at end of file diff --git a/duneana/MCParticleAna/fcl/mcparticle_ana_gamma_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/MCParticleAna/fcl/mcparticle_ana_gamma_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..780fde17 --- /dev/null +++ b/duneana/MCParticleAna/fcl/mcparticle_ana_gamma_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,4 @@ +#include "mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.mcparticleana.MCParticlePDG: 22 # For gamma particles +physics.analyzers.mcparticleana.MCParticleMinKE: 4. # MeV \ No newline at end of file diff --git a/duneana/MCParticleAna/fcl/mcparticle_ana_neutron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/MCParticleAna/fcl/mcparticle_ana_neutron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..dfeac02b --- /dev/null +++ b/duneana/MCParticleAna/fcl/mcparticle_ana_neutron_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,3 @@ +#include "mcparticle_ana_all_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.mcparticleana.MCParticlePDG: 2112 # For neutron particles \ No newline at end of file diff --git a/duneana/SolarNuAna/SolarNuAna_module.cc b/duneana/SolarNuAna/SolarNuAna_module.cc index 85e07ec1..01bf6ca1 100644 --- a/duneana/SolarNuAna/SolarNuAna_module.cc +++ b/duneana/SolarNuAna/SolarNuAna_module.cc @@ -79,22 +79,25 @@ 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; - std::string fSignalLabel, fClusterLabel, fSolarClusterLabel, fClusterChargeVariable, fOpHitTimeVariable, fAdjOpFlashMinPEAttenuate; - int fClusterAlgoAdjChannel, fClusterInd0MatchTime, fClusterInd1MatchTime, fClusterPreselectionNHits, fAdjOpFlashMinNHitCut, fAdjOpFlashMinPEAttenuationStrength; + std::string fSignalLabel, fClusterLabel, fSolarClusterLabel, fClusterChargeVariable, fOpHitTimeVariable, fAdjOpFlashMinPEAttenuate, fAdjOpFlashMaxPEAttenuate, fFlashMatchBy; + int fClusterAlgoAdjChannel, fClusterInd0MatchTime, fClusterInd1MatchTime, fClusterPreselectionNHits, fAdjOpFlashMinNHitCut, fAdjOpFlashMinPEAttenuationStrength, fAdjOpFlashMaxPEAttenuationStrength; 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; + float fClusterMatchTime, fAdjClusterRad, fMinClusterCharge, fClusterMatchCharge, fClusterMatchNHit, fClusterAlgoTime; + float fOpFlashTimeOffset, fOpFlashAlgoMinTime, fOpFlashAlgoMaxTime, fOpFlashAlgoRad, fOpFlashAlgoPE, fOpFlashAlgoTriggerPE, fOpFlashAlgoHotVertexThld; + float fAdjOpFlashX, fAdjOpFlashY, fAdjOpFlashZ, fAdjOpFlashMaxPERatioCut, fAdjOpFlashMinPECut, fAdjOpFlashMaxPECut, fAdjOpFlashMinPEAttenuation, fAdjOpFlashMaxPEAttenuation; + float fXACathodeX, fXAMembraneY, fXAStartCapZ, fXAFinalCapZ, fFlashMatchByPELightMapExponent; + bool fOpFlashAlgoHitDuplicates, fOpFlashAlgoWeightedTime; bool fClusterPreselectionSignal, fClusterPreselectionPrimary, fClusterPreselectionTrack, fClusterPreselectionFlashMatch; - bool fGenerateSolarCluster, fGenerateAdjCluster, fGenerateAdjOpFlash, fFlashMatchByResidual; + bool fGenerateSolarCluster, fGenerateAdjCluster, fGenerateAdjOpFlash; 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. + int AnalyzedEvents = 0; // Total number of analyzed events. // --- Our TTrees, and its associated variables. TTree *fConfigTree; TTree *fMCTruthTree; @@ -119,12 +122,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 @@ -148,7 +154,7 @@ namespace solar std::unique_ptr adjophits; std::unique_ptr lowe; }; -#endif +#endif // SolarNuAna_h //...................................................... SolarNuAna::SolarNuAna(fhicl::ParameterSet const &p) @@ -171,13 +177,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 +206,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,16 +216,21 @@ 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); + fAdjOpFlashMinNHitCut = p.get("AdjOpFlashMinNHitCut"); + 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"); + fAdjOpFlashMinPECut = p.get("AdjOpFlashMinPECut", 0.0); + fAdjOpFlashMaxPECut = p.get("AdjOpFlashMaxPECut", 1e9); fAdjOpFlashMinPEAttenuate = p.get("AdjOpFlashMinPEAttenuate"); - fAdjOpFlashMinPEAttenuation = p.get("AdjOpFlashMinPEAttenuation", 0.9); - fAdjOpFlashMinPEAttenuationStrength = p.get("AdjOpFlashMinPEAttenuationStrength", 4); - fAdjOpFlashMinNHitCut = p.get("AdjOpFlashMinNHitCut"); - fFlashMatchByResidual = p.get("FlashMatchByResidual"); + fAdjOpFlashMaxPEAttenuate = p.get("AdjOpFlashMaxPEAttenuate"); + fAdjOpFlashMinPEAttenuation = p.get("AdjOpFlashMinPEAttenuation", 0.0); + fAdjOpFlashMaxPEAttenuation = p.get("AdjOpFlashMaxPEAttenuation", 0.0); + fAdjOpFlashMinPEAttenuationStrength = p.get("AdjOpFlashMinPEAttenuationStrength", 0); + fAdjOpFlashMaxPEAttenuationStrength = p.get("AdjOpFlashMaxPEAttenuationStrength", 0); + fFlashMatchBy = p.get("FlashMatchBy", "maximum"); + fFlashMatchByPELightMapExponent = p.get("FlashMatchByPELightMapExponent", 1); fSaveSignalDaughters = p.get("SaveSignalDaughters"); fSaveSignalEDep = p.get("SaveSignalEDep"); fSaveSignalOpHits = p.get("SaveSignalOpHits"); @@ -250,6 +264,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 +289,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); @@ -281,16 +297,21 @@ namespace solar fConfigTree->Branch("OpFlashAlgoHitDuplicates", &fOpFlashAlgoHitDuplicates); fConfigTree->Branch("AdjOpFlashMembraneProjection", &fAdjOpFlashMembraneProjection); fConfigTree->Branch("AdjOpFlashEndCapProjection", &fAdjOpFlashEndCapProjection); + fConfigTree->Branch("AdjOpFlashMinNHitCut", &fAdjOpFlashMinNHitCut); fConfigTree->Branch("AdjOpFlashX", &fAdjOpFlashX); fConfigTree->Branch("AdjOpFlashY", &fAdjOpFlashY); fConfigTree->Branch("AdjOpFlashZ", &fAdjOpFlashZ); fConfigTree->Branch("AdjOpFlashMaxPERatioCut", &fAdjOpFlashMaxPERatioCut); fConfigTree->Branch("AdjOpFlashMinPECut", &fAdjOpFlashMinPECut); + fConfigTree->Branch("AdjOpFlashMaxPECut", &fAdjOpFlashMaxPECut); fConfigTree->Branch("AdjOpFlashMinPEAttenuate", &fAdjOpFlashMinPEAttenuate); + fConfigTree->Branch("AdjOpFlashMaxPEAttenuate", &fAdjOpFlashMaxPEAttenuate); fConfigTree->Branch("AdjOpFlashMinPEAttenuation", &fAdjOpFlashMinPEAttenuation); + fConfigTree->Branch("AdjOpFlashMaxPEAttenuation", &fAdjOpFlashMaxPEAttenuation); fConfigTree->Branch("AdjOpFlashMinPEAttenuationStrength", &fAdjOpFlashMinPEAttenuationStrength); - fConfigTree->Branch("AdjOpFlashMinNHitCut", &fAdjOpFlashMinNHitCut); - fConfigTree->Branch("FlashMatchByResidual", &fFlashMatchByResidual); + fConfigTree->Branch("AdjOpFlashMaxPEAttenuationStrength", &fAdjOpFlashMaxPEAttenuationStrength); + fConfigTree->Branch("FlashMatchBy", &fFlashMatchBy); + fConfigTree->Branch("FlashMatchByPELightMapExponent", &fFlashMatchByPELightMapExponent); fConfigTree->Branch("SaveSignalDaughters", &fSaveSignalDaughters); fConfigTree->Branch("SaveSignalEDep", &fSaveSignalEDep); fConfigTree->Branch("SaveSignalOpHits", &fSaveSignalOpHits); @@ -298,6 +319,7 @@ namespace solar fConfigTree->Branch("SaveAdjOpFlashInfo", &fSaveAdjOpFlashInfo); fConfigTree->Branch("SaveTrackInfo", &fSaveTrackInfo); fConfigTree->Branch("SelectedEvents", &SelectedEvents); + fConfigTree->Branch("AnalyzedEvents", &AnalyzedEvents); // MC Truth info. fMCTruthTree->Branch("Event", &Event, "Event/I"); // Event number @@ -311,7 +333,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 +346,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 +392,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 +408,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 +435,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 +458,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 +485,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 +511,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); @@ -526,7 +554,7 @@ namespace solar hDriftTime = tfs->make("hDriftTime", "hDriftTime", 100, -400., 400., 100, 0., 10000.); hXTruth = tfs->make("hXTruth", "Missmatch in X distance; Distance [cm]; True X position [cm]", 100, -600, 600, 100, -600, 600); hYTruth = tfs->make("hYTruth", "Missmatch in Y distance; Distance [cm]; True Y position [cm]", 100, -600, 600, 100, -600, 600); - hZTruth = tfs->make("hZTruth", "Missmatch in Z distance; Distance [cm]; True Z position [cm]", 100, -600, 600, 100, 0, 1600); + hZTruth = tfs->make("hZTruth", "Missmatch in Z distance; Distance [cm]; True Z position [cm]", 100, -600, 600, 100, 0, 2100); } // BeginJob //...................................................... @@ -732,7 +760,7 @@ namespace solar } sSignalTruth = sSignalTruth + sSignalParticle + " Energy: " + ProducerUtils::str(SignalParticleK) + " MeV\n"; - sSignalTruth += "\t- Strat Position: (" + ProducerUtils::str(SignalParticleX) + ", " + ProducerUtils::str(SignalParticleY) + ", " + ProducerUtils::str(SignalParticleZ) + ") cm\n"; + sSignalTruth += "\t- Start Position: (" + ProducerUtils::str(SignalParticleX) + ", " + ProducerUtils::str(SignalParticleY) + ", " + ProducerUtils::str(SignalParticleZ) + ") cm\n"; sSignalTruth += "\t- Final Position: (" + ProducerUtils::str(SignalParticle.EndX()) + ", " + ProducerUtils::str(SignalParticle.EndY()) + ", " + ProducerUtils::str(SignalParticle.EndZ()) + ") cm\n"; } } @@ -880,954 +908,1011 @@ namespace solar producer->PrintInColor("\nKinetic energy of signal particle is above threshold of " + ProducerUtils::str(fMaxSignalK) + " MeV. Skipping event.\n", ProducerUtils::GetColor("red"), "Warning"); } if (SelectedEvent) { - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - //---------------------------------------------------------------------- PMTrack Analysis -----------------------------------------------------------------------// - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - art::Handle> TrackHandle; - std::vector> TrackList; - if (evt.getByLabel(fTrackLabel, TrackHandle)) { - art::fill_ptr_vector(TrackList, TrackHandle); - } - TrackNum = int(TrackList.size()); - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - //------------------------------------------------------------------- Optical Flash Analysis --------------------------------------------------------------------// - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - // Find OpHits and OpFlashes associated with the event - std::string sOpFlashTruth = ""; - std::vector> OpHitList; - art::Handle> OpHitHandle; - std::vector>> OpHitVec; - std::vector> OpHitIdx; - if (evt.getByLabel(fOpHitLabel, OpHitHandle)) { - art::fill_ptr_vector(OpHitList, OpHitHandle); - } - // Grab assns with OpHits to get match to neutrino purity - OpHitNum = int(OpHitList.size()); - if (fGenerateAdjOpFlash) { - fOpFlashLabel = "solarflash"; - std::vector FlashVec; - 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; - OpFlashPlane.push_back(TheFlash.Plane); - OpFlashNHits.push_back(TheFlash.NHit); - OpFlashTime.push_back(TheFlash.Time - fOpFlashTimeOffset); // Convert to microseconds happens in AdjOpHits - OpFlashDeltaT.push_back(TheFlash.TimeWidth); // Convert to microseconds - OpFlashPE.push_back(TheFlash.PE); - OpFlashMaxPE.push_back(TheFlash.MaxPE); - 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); - 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; - for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) - { - if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) - ThisOphitPurity += 1; - } - // Check if ThisOpHitTrackIds is empty - if (ThisOpHitTrackIds.size() == 0) - ThisOphitPurity = 0; - else - ThisOphitPurity /= int(ThisOpHitTrackIds.size()); - - ThisOpFlashPur += ThisOphitPurity * OpHit.PE(); - auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(OpHit.OpChannel()).GetCenter(); - SOpHitPur.push_back(ThisOphitPurity); - SOpHitChannel.push_back(OpHit.OpChannel()); - - if (fOpHitTimeVariable == "StartTime") - 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 - - 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); - } - // Check if OpHitVec[i] is empty - if (OpHitVec[i].size() == 0) - ThisOpFlashPur = 0; - else - ThisOpFlashPur /= TheFlash.PE; - - OpFlashPur.push_back(ThisOpFlashPur); - 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 += " - 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*** 2nd Sanity check: #OpHits " + ProducerUtils::str(int(OpHitVec[i].size())) + " >= " + ProducerUtils::str(TheFlash.NHit) + "\n"; - } + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + //---------------------------------------------------------------------- PMTrack Analysis -----------------------------------------------------------------------// + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + art::Handle> TrackHandle; + std::vector> TrackList; + if (evt.getByLabel(fTrackLabel, TrackHandle)) { + art::fill_ptr_vector(TrackList, TrackHandle); } - } - else { - std::vector> OpFlashList; - art::Handle> FlashHandle; - - if (evt.getByLabel(fOpFlashLabel, FlashHandle)){ - art::fill_ptr_vector(OpFlashList, FlashHandle); + TrackNum = int(TrackList.size()); + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + //------------------------------------------------------------------- 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; + std::vector>> OpHitVec; + std::vector> OpHitIdx; + if (evt.getByLabel(fOpHitLabel, OpHitHandle)) { + art::fill_ptr_vector(OpHitList, OpHitHandle); } // Grab assns with OpHits to get match to neutrino purity - OpFlashNum = int(OpFlashList.size()); - art::FindManyP OpAssns(OpFlashList, evt, fOpFlashLabel); - // Loop over OpFlashList and assign OpHits to each flash - for (int i = 0; i < int(OpFlashList.size()); i++) - { - recob::OpFlash TheFlash = *OpFlashList[i]; - std::vector> MatchedHits = OpAssns.at(i); - int NMatchedHits = MatchedHits.size(); - double FlashStdDev = 0.0, TotalFlashPE = 0, MaxOpHitPE = 0; - std::vector varXY, varYZ, varXZ; - varXY = varYZ = varXZ = {}; - - for (int j = 0; j < NMatchedHits; j++) - { // Loop over OpHits in the flash - recob::OpHit OpHit = *MatchedHits[j]; - art::Ptr OpHitPtr = MatchedHits[j]; - mf::LogDebug("SolarNuAna") << "Assigning OpHit to Flash"; - const std::vector ThisOpHitTrackIds = pbt->OpHitToTrackIds(OpHit); - float ThisOphitPurity = 0; + OpHitNum = int(OpHitList.size()); + if (fGenerateAdjOpFlash) { + fOpFlashLabel = "solarflash"; + std::vector FlashVec; + 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; - for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) - { - if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) { - ThisOphitPurity += 1; - } - } - - auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(OpHit.OpChannel()).GetCenter(); - TotalFlashPE += OpHit.PE(); + if (fOpFlashAlgoWeightedTime) + ThisOpFlashTime = TheFlash.TimeWeighted - fOpFlashTimeOffset; + else + ThisOpFlashTime = TheFlash.MainOpHitTime - fOpFlashTimeOffset; - 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())); + OpHitAmplitude.push_back(TheFlash.MainOpHitAmplitude); + OpFlashPlane.push_back(TheFlash.Plane); + OpFlashNHits.push_back(TheFlash.NHit); + OpFlashTime.push_back(ThisOpFlashTime); + OpFlashDeltaT.push_back(TheFlash.TimeWidth); // Convert to microseconds + OpFlashPE.push_back(TheFlash.PE); + 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); - if (OpHit.PE() > MaxOpHitPE) { - MaxOpHitPE = OpHit.PE(); - } + for (int j = 0; j < int(OpHitVec[i].size()); j++) + { + recob::OpHit OpHit = *OpHitVec[i][j]; - SOpHitFlashID.push_back(i); - SOpHitPE.push_back(OpHit.PE()); - SOpHitX.push_back(OpHitXYZ.X()); - SOpHitY.push_back(OpHitXYZ.Y()); - SOpHitZ.push_back(OpHitXYZ.Z()); - - if (fOpHitTimeVariable == "StartTime") - 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 + const std::vector ThisOpHitTrackIds = pbt->OpHitToTrackIds(OpHit); + float ThisOpHitTime = -1e6; + float ThisOpHitPurity = 0; + for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) + { + if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) + ThisOpHitPurity += 1; + } + // Check if ThisOpHitTrackIds is empty + if (ThisOpHitTrackIds.size() == 0) + ThisOpHitPurity = 0; + else + ThisOpHitPurity /= int(ThisOpHitTrackIds.size()); + + ThisOpFlashPur += ThisOpHitPurity * OpHit.PE(); + auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(OpHit.OpChannel()).GetCenter(); + SOpHitPur.push_back(ThisOpHitPurity); + SOpHitChannel.push_back(OpHit.OpChannel()); + + if (fOpHitTimeVariable == "StartTime") + ThisOpHitTime = OpHit.StartTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + else // Default to PeakTime + 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 - SOpHitChannel.push_back(OpHit.OpChannel()); - SOpHitPlane.push_back(adjophits->GetOpHitPlane(OpHitPtr, 0.1)); // Get plane assignment for the OpHit - } // End of OpHit loop - - OpHitVec.push_back(MatchedHits); - FlashStdDev = adjophits->GetOpFlashPlaneSTD(TheFlash.Frame(), varXY, varYZ, varXZ); - int TerminalOutput = ProducerUtils::supress_stdout(); - double ThisOpFlashPur = pbt->OpHitCollectionPurity(SignalTrackIDs, MatchedHits); - ProducerUtils::resume_stdout(TerminalOutput); - - // Calculate the flash purity, only for the Signal events - OpFlashID.push_back(i); - OpFlashPlane.push_back(TheFlash.Frame()); - OpFlashPur.push_back(ThisOpFlashPur); - OpFlashMaxPE.push_back(MaxOpHitPE); - OpFlashSTD.push_back(FlashStdDev); - OpFlashX.push_back(TheFlash.XCenter()); - OpFlashY.push_back(TheFlash.YCenter()); - OpFlashZ.push_back(TheFlash.ZCenter()); - OpFlashPE.push_back(TheFlash.TotalPE()); - OpFlashNHits.push_back(MatchedHits.size()); - OpFlashFast.push_back(TheFlash.FastToTotal()); - - 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 + 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) + ThisOpFlashPur = 0; + else + 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(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.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 { - OpFlashTime.push_back(TheFlash.Time()); // Expected flash to provide time in microseconds - OpFlashDeltaT.push_back(TheFlash.TimeWidth()); // Expected flash to provide time width in microseconds + } + 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; + + if (evt.getByLabel(fOpFlashLabel, FlashHandle)){ + art::fill_ptr_vector(OpFlashList, FlashHandle); } + // Grab assns with OpHits to get match to neutrino purity + OpFlashNum = int(OpFlashList.size()); + art::FindManyP OpAssns(OpFlashList, evt, fOpFlashLabel); + // Loop over OpFlashList and assign OpHits to each flash + 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; + varXY = varYZ = varXZ = {}; + + for (int j = 0; j < NMatchedHits; j++) + { // Loop over OpHits in the flash + recob::OpHit OpHit = *MatchedHits[j]; + art::Ptr OpHitPtr = MatchedHits[j]; + mf::LogDebug("SolarNuAna") << "Assigning OpHit to Flash"; + const std::vector ThisOpHitTrackIds = pbt->OpHitToTrackIds(OpHit); + float ThisOpHitPurity = 0; + float ThisOpHitTime = -1e6; + + for (auto const &ThisOpHitTrackId : ThisOpHitTrackIds) + { + if (SignalTrackIDs.find(ThisOpHitTrackId) != SignalTrackIDs.end()) { + ThisOpHitPurity += 1; + } + } - 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"; - sOpFlashTruth += " - Vertex (" + ProducerUtils::str(TheFlash.XCenter()) + ", " + ProducerUtils::str(TheFlash.YCenter()) + ", " + ProducerUtils::str(TheFlash.ZCenter()) + ")\n"; - sOpFlashTruth += "\t*** 1st Sanity check: Ratio " + ProducerUtils::str(MaxOpHitPE / TotalFlashPE) + " <= " + ProducerUtils::str(fAdjOpFlashMaxPERatioCut) + "\n"; - sOpFlashTruth += "\t*** 2nd Sanity check: #OpHits " + ProducerUtils::str(int(NMatchedHits)) + " >= " + ProducerUtils::str(int(TheFlash.PEs().size())) + "\n"; - } - } - } - sOpFlashTruth = sOpFlashTruth + "\n# of OpFlashes (" + fOpFlashLabel + ") in full geometry: " + ProducerUtils::str(OpFlashNum) + "\n"; - producer->PrintInColor(sOpFlashTruth, ProducerUtils::GetColor("blue")); + auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(OpHit.OpChannel()).GetCenter(); + TotalFlashPE += OpHit.PE(); + + 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())); + + 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; + } + } - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - //---------------------------------------------------------------- Hit collection and assignment ----------------------------------------------------------------// - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - // --- Lift out the reco hits: - auto RecoHits = evt.getValidHandle>(fHitLabel); - std::vector> RecoHitsPtr; - int NTotHits = RecoHits->size(); + SOpHitFlashID.push_back(i); + SOpHitPE.push_back(OpHit.PE()); + SOpHitX.push_back(OpHitXYZ.X()); + SOpHitY.push_back(OpHitXYZ.Y()); + SOpHitZ.push_back(OpHitXYZ.Z()); + + if (fOpHitTimeVariable == "StartTime") + ThisOpHitTime = OpHit.StartTime(); // Convert to microseconds + // SOpHitTime.push_back(OpHit.StartTime() * clockData.OpticalClock().TickPeriod() - fOpFlashTimeOffset); // Convert to microseconds + else // Default to PeakTime + 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); + FlashStdDev = adjophits->GetOpFlashPlaneSTD(TheFlash.Frame(), varXY, varYZ, varXZ); + int TerminalOutput = ProducerUtils::supress_stdout(); + double ThisOpFlashPur = pbt->OpHitCollectionPurity(SignalTrackIDs, MatchedHits); + ProducerUtils::resume_stdout(TerminalOutput); - for (int i = 0; i < NTotHits; ++i) - { - // --- Loop over the reconstructed hits to separate them among tpc planes according to view and signal type - recob::Hit const &ThisHit = RecoHits->at(i); - // Add to RecoHitsPtr - RecoHitsPtr.push_back(art::Ptr(RecoHits, i)); - if (ThisHit.View() == 0) - { - Ind0Hits.push_back(ThisHit); - } // SignalType = 0 - else if (ThisHit.View() == 1) - { - Ind1Hits.push_back(ThisHit); - } // SignalType = 0 - else if (ThisHit.View() == 2) - { - ColHits.push_back(ThisHit); - } // SignalType = 1 - else - { - GhostHits.push_back(ThisHit); - mf::LogError("SolarNuAna") << "Hit was found with view out of scope"; - } - } + // Calculate the flash purity, only for the Signal events + OpFlashID.push_back(i); + OpFlashPlane.push_back(TheFlash.Frame()); + OpFlashPur.push_back(ThisOpFlashPur); + OpFlashMaxPE.push_back(MaxOpHitPE); + OpFlashSTD.push_back(FlashStdDev); + OpFlashX.push_back(TheFlash.XCenter()); + OpFlashY.push_back(TheFlash.YCenter()); + OpFlashZ.push_back(TheFlash.ZCenter()); + 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() - fOpFlashTimeOffset); // Expected flash to provide time in microseconds + OpFlashDeltaT.push_back(TheFlash.TimeWidth()); // Expected flash to provide time width in microseconds + } - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - //-------------------------------------------------------------- Cluster creation and analysis ------------------------------------------------------------------// - //---------------------------------------------------------------------------------------------------------------------------------------------------------------// - std::string sRecoObjects = ""; - std::vector>> ClustersPtr; - std::vector PerPlaneClusters; - std::vector>> AllPlaneClusters; - std::vector> ClustersIdx = {{}, {}, {}}; - std::vector> RecoHitIdx; - // Map to associate the ClusterIdx with the position in the ClVectors - std::map> ClIdxMap; - - if (fGenerateSolarCluster == false) { - // Get clusters from event recob::Cluster with label "planecluster" - // ... - } - else { - lowe->CalcAdjHits(RecoHitsPtr, ClustersPtr, RecoHitIdx, evt); - for (int i = 0; i < int(ClustersPtr.size()); i++) + 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"; + sOpFlashTruth += " - Vertex (" + ProducerUtils::str(TheFlash.XCenter()) + ", " + ProducerUtils::str(TheFlash.YCenter()) + ", " + ProducerUtils::str(TheFlash.ZCenter()) + ")\n"; + sOpFlashTruth += "\t*** 1st Sanity check: Ratio " + ProducerUtils::str(MaxOpHitPE / TotalFlashPE) + " <= " + ProducerUtils::str(fAdjOpFlashMaxPERatioCut) + "\n"; + sOpFlashTruth += "\t*** 2nd Sanity check: #OpHits " + ProducerUtils::str(int(NMatchedHits)) + " >= " + ProducerUtils::str(int(TheFlash.PEs().size())) + "\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")); + + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + //---------------------------------------------------------------- Hit collection and assignment ----------------------------------------------------------------// + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + // --- Lift out the reco hits: + auto RecoHits = evt.getValidHandle>(fHitLabel); + std::vector> RecoHitsPtr; + int NTotHits = RecoHits->size(); + + for (int i = 0; i < NTotHits; ++i) { - std::vector ThisHitVector = {}; // Convert pointer to vector - for (int j = 0; j < int(ClustersPtr[i].size()); j++) + // --- Loop over the reconstructed hits to separate them among tpc planes according to view and signal type + recob::Hit const &ThisHit = RecoHits->at(i); + // Add to RecoHitsPtr + RecoHitsPtr.push_back(art::Ptr(RecoHits, i)); + if (ThisHit.View() == 0) { - ThisHitVector.push_back(*ClustersPtr[i][j]); - } - int ThisIdx = RecoHitIdx[i][0]; - if (RecoHitsPtr[ThisIdx]->View() == 0) + Ind0Hits.push_back(ThisHit); + } // SignalType = 0 + else if (ThisHit.View() == 1) { - Clusters0.push_back(ThisHitVector); - ClustersIdx[0].push_back(i); - } - else if (RecoHitsPtr[ThisIdx]->View() == 1) + Ind1Hits.push_back(ThisHit); + } // SignalType = 0 + else if (ThisHit.View() == 2) { - Clusters1.push_back(ThisHitVector); - ClustersIdx[1].push_back(i); - } - else if (RecoHitsPtr[ThisIdx]->View() == 2) + ColHits.push_back(ThisHit); + } // SignalType = 1 + else { - Clusters2.push_back(ThisHitVector); - ClustersIdx[2].push_back(i); + GhostHits.push_back(ThisHit); + mf::LogError("SolarNuAna") << "Hit was found with view out of scope"; } - else if (RecoHitsPtr[ThisIdx]->View() == 3) + } + + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + //-------------------------------------------------------------- Cluster creation and analysis ------------------------------------------------------------------// + //---------------------------------------------------------------------------------------------------------------------------------------------------------------// + std::string sRecoObjects = ""; + std::vector>> ClustersPtr; + std::vector PerPlaneClusters; + std::vector>> AllPlaneClusters; + std::vector> ClustersIdx = {{}, {}, {}}; + std::vector> RecoHitIdx; + // Map to associate the ClusterIdx with the position in the ClVectors + std::map> ClIdxMap; + + if (fGenerateSolarCluster == false) { + // Get clusters from event recob::Cluster with label "planecluster" + // ... + } + else { + lowe->CalcAdjHits(RecoHitsPtr, ClustersPtr, RecoHitIdx, evt); + for (int i = 0; i < int(ClustersPtr.size()); i++) { - Clusters3.push_back(ThisHitVector); + std::vector ThisHitVector = {}; // Convert pointer to vector + for (int j = 0; j < int(ClustersPtr[i].size()); j++) + { + ThisHitVector.push_back(*ClustersPtr[i][j]); + } + int ThisIdx = RecoHitIdx[i][0]; + if (RecoHitsPtr[ThisIdx]->View() == 0) + { + Clusters0.push_back(ThisHitVector); + ClustersIdx[0].push_back(i); + } + else if (RecoHitsPtr[ThisIdx]->View() == 1) + { + Clusters1.push_back(ThisHitVector); + ClustersIdx[1].push_back(i); + } + else if (RecoHitsPtr[ThisIdx]->View() == 2) + { + Clusters2.push_back(ThisHitVector); + ClustersIdx[2].push_back(i); + } + else if (RecoHitsPtr[ThisIdx]->View() == 3) + { + Clusters3.push_back(ThisHitVector); + } } + lowe->MakeClusterVector(PerPlaneClusters, ClustersPtr, evt); + AllPlaneClusters = {Clusters0, Clusters1, Clusters2}; } - lowe->MakeClusterVector(PerPlaneClusters, ClustersPtr, evt); - AllPlaneClusters = {Clusters0, Clusters1, Clusters2}; - } - - std::vector>> ClVecGenPur = {{}, {}, {}}; - std::vector> ClMainID = {{}, {}, {}}, ClTPC = {{}, {}, {}}, ClNHits = {{}, {}, {}}, ClGen = {{}, {}, {}}; - std::vector> ClCharge = {{}, {}, {}}, ClMaxCharge = {{}, {}, {}}, ClT = {{}, {}, {}}, ClX = {{}, {}, {}}, ClY = {{}, {}, {}}, ClZ = {{}, {}, {}}; - std::vector> ClFracE = {{}, {}, {}}, ClFracGa = {{}, {}, {}}, ClFracNe = {{}, {}, {}}, ClFracRest = {{}, {}, {}}; - std::vector> ClPur = {{}, {}, {}}, Cldzdy = {{}, {}, {}}, ClGenPur = {{}, {}, {}}; - - sRecoObjects += "\n# OpHits (" + fOpHitLabel + ") in full geometry: " + ProducerUtils::str(OpHitNum); - sRecoObjects += "\n# OpFlashes (" + fOpFlashLabel + ") in full geometry: " + ProducerUtils::str(OpFlashNum); - sRecoObjects += "\n# Hits (" + fHitLabel + ") in each view: " + ProducerUtils::str(int(Ind0Hits.size())) + ", " + ProducerUtils::str(int(Ind1Hits.size())) + ", " + ProducerUtils::str(int(ColHits.size())) + ", " + ProducerUtils::str(int(GhostHits.size())); - sRecoObjects += "\n# Cluster from the hits: " + ProducerUtils::str(int(Clusters0.size())) + ", " + ProducerUtils::str(int(Clusters1.size())) + ", " + ProducerUtils::str(int(Clusters2.size())) + ", " + ProducerUtils::str(int(Clusters3.size())); - sRecoObjects += "\n# Tracks (" + fTrackLabel + ") in full geometry: " + ProducerUtils::str(TrackNum); - producer->PrintInColor(sRecoObjects, ProducerUtils::GetColor("cyan")); - - //------------------------------------------------------------ First complete cluster analysis ------------------------------------------------------------------// - // --- Now loop over the planes and the clusters to calculate the cluster properties - for (int idx = 0; idx < 3; idx++) - { - int nhit, clustTPC; - float FracE, FracGa, FracNe, FracRest, clustX, clustY, clustZ, clustT, ncharge, maxHit, dzdy; - std::vector> Clusters = AllPlaneClusters[idx]; - // --- Loop over the clusters - for (int i = 0; i < int(Clusters.size()); i++) + 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 = {{}, {}, {}}; + std::vector> ClFracE = {{}, {}, {}}, ClFracGa = {{}, {}, {}}, ClFracNe = {{}, {}, {}}, ClFracRest = {{}, {}, {}}; + std::vector> ClPur = {{}, {}, {}}, Cldzdy = {{}, {}, {}}, ClGenPur = {{}, {}, {}}; + + sRecoObjects += "\n# OpHits (" + fOpHitLabel + ") in full geometry: " + ProducerUtils::str(OpHitNum); + sRecoObjects += "\n# OpFlashes (" + fOpFlashLabel + ") in full geometry: " + ProducerUtils::str(OpFlashNum); + sRecoObjects += "\n# Hits (" + fHitLabel + ") in each view: " + ProducerUtils::str(int(Ind0Hits.size())) + ", " + ProducerUtils::str(int(Ind1Hits.size())) + ", " + ProducerUtils::str(int(ColHits.size())) + ", " + ProducerUtils::str(int(GhostHits.size())); + sRecoObjects += "\n# Cluster from the hits: " + ProducerUtils::str(int(Clusters0.size())) + ", " + ProducerUtils::str(int(Clusters1.size())) + ", " + ProducerUtils::str(int(Clusters2.size())) + ", " + ProducerUtils::str(int(Clusters3.size())); + sRecoObjects += "\n# Tracks (" + fTrackLabel + ") in full geometry: " + ProducerUtils::str(TrackNum); + producer->PrintInColor(sRecoObjects, ProducerUtils::GetColor("cyan")); + + //------------------------------------------------------------ First complete cluster analysis ------------------------------------------------------------------// + // --- Now loop over the planes and the clusters to calculate the cluster properties + for (int idx = 0; idx < 3; idx++) { - int MainTrID; - int MainGenerator = 0; - float Pur = 0; - std::vector thisdzdy = {}; - - nhit = Clusters[i].size(); - ncharge = maxHit = clustT = FracE = FracGa = FracNe = FracRest = clustX = clustY = clustZ = dzdy = 0; - clustTPC = -1; // Initialize clustTPC to -1 - // Define a vector of floats with size equal to the number of generators + 1 - std::vector VecGenPur(fLabels.size() + 1, 0); - - for (recob::Hit TPCHit : Clusters[i]) + int nhit, clustTPC; + float FracE, FracGa, FracNe, FracRest, clustX, clustY, clustZ, clustT, ncharge, maxHit, dzdy; + std::vector> Clusters = AllPlaneClusters[idx]; + + // --- Loop over the clusters + for (int i = 0; i < int(Clusters.size()); i++) { - ncharge += TPCHit.Integral(); - const geo::WireGeo *wire = wireReadout.WirePtr(TPCHit.WireID()); // Wire directions should be the same for all hits of the same view (can be used to check) - double hitCharge; - - geo::Point_t hXYZ = wire->GetCenter(); - geo::Point_t sXYZ = wire->GetStart(); - geo::Point_t eXYZ = wire->GetEnd(); - geo::Vector_t direction = eXYZ - sXYZ; - auto dyds = direction.Y(), dzds = direction.Z(); - thisdzdy.push_back(dzds / dyds); - - int TPC = TPCHit.WireID().TPC; - clustX += TPCHit.Integral() * hXYZ.X(); - clustY += TPCHit.Integral() * hXYZ.Y(); - clustZ += TPCHit.Integral() * hXYZ.Z(); - clustT += TPCHit.Integral() * TPCHit.PeakTime() * clockData.TPCClock().TickPeriod(); // Convert to microseconds - - if (TPCHit.Integral() > maxHit) { // If clusterTPC not in TPCIDMap, set it to -1 - if (TPCIDMap.find(TPC) == TPCIDMap.end()) { - clustTPC = -1; + int MainTrID; + int MainGenerator = 0; + float Pur = 0; + std::vector thisdzdy = {}; + + nhit = Clusters[i].size(); + ncharge = maxHit = clustT = FracE = FracGa = FracNe = FracRest = clustX = clustY = clustZ = dzdy = 0; + clustTPC = -1; // Initialize clustTPC to -1 + // Define a vector of floats with size equal to the number of generators + 1 + std::vector VecGenPur(fLabels.size() + 1, 0); + + for (recob::Hit TPCHit : Clusters[i]) + { + ncharge += TPCHit.Integral(); + const geo::WireGeo *wire = wireReadout.WirePtr(TPCHit.WireID()); // Wire directions should be the same for all hits of the same view (can be used to check) + double hitCharge; + + geo::Point_t hXYZ = wire->GetCenter(); + geo::Point_t sXYZ = wire->GetStart(); + geo::Point_t eXYZ = wire->GetEnd(); + geo::Vector_t direction = eXYZ - sXYZ; + auto dyds = direction.Y(), dzds = direction.Z(); + thisdzdy.push_back(dzds / dyds); + + int TPC = TPCHit.WireID().TPC; + clustX += TPCHit.Integral() * hXYZ.X(); + clustY += TPCHit.Integral() * hXYZ.Y(); + clustZ += TPCHit.Integral() * hXYZ.Z(); + clustT += TPCHit.Integral() * TPCHit.PeakTime() * clockData.TPCClock().TickPeriod(); // Convert to microseconds + + if (TPCHit.Integral() > maxHit) { // If clusterTPC not in TPCIDMap, set it to -1 + if (TPCIDMap.find(TPC) == TPCIDMap.end()) { + clustTPC = -1; + } + else { + clustTPC = TPC; + } + // Look for maxHit inside cluster + maxHit = TPCHit.Integral(); } else { clustTPC = TPC; } - // Look for maxHit inside cluster - maxHit = TPCHit.Integral(); - } - else { - clustTPC = TPC; - } - MainTrID = 0; - double TopEFrac = 0; - std::vector ThisHitIDE = bt_serv->HitToTrackIDEs(clockData, TPCHit); + MainTrID = 0; + double TopEFrac = 0; + std::vector ThisHitIDE = bt_serv->HitToTrackIDEs(clockData, TPCHit); - for (size_t ideL = 0; ideL < ThisHitIDE.size(); ++ideL) - { - if (ThisHitIDE[ideL].energyFrac > TopEFrac) { - TopEFrac = ThisHitIDE[ideL].energyFrac; - MainTrID = abs(ThisHitIDE[ideL].trackID); + for (size_t ideL = 0; ideL < ThisHitIDE.size(); ++ideL) + { + if (ThisHitIDE[ideL].energyFrac > TopEFrac) { + TopEFrac = ThisHitIDE[ideL].energyFrac; + MainTrID = abs(ThisHitIDE[ideL].trackID); + } } - } - for (int frac = 0; frac < int(ClPartTrackIDs.size()); ++frac) - { - for (int trck = 0; trck < int(ClPartTrackIDs[frac].size()); ++trck) + for (int frac = 0; frac < int(ClPartTrackIDs.size()); ++frac) { - if (MainTrID == ClPartTrackIDs[frac][trck]) { - if (frac == 0) { - FracE = FracE + TPCHit.Integral(); - } - else if (frac == 1) { - FracGa = FracGa + TPCHit.Integral(); - } - else if (frac == 2) { - FracNe = FracNe + TPCHit.Integral(); - } - else { - FracRest = FracRest + TPCHit.Integral(); + for (int trck = 0; trck < int(ClPartTrackIDs[frac].size()); ++trck) + { + if (MainTrID == ClPartTrackIDs[frac][trck]) { + if (frac == 0) { + FracE = FracE + TPCHit.Integral(); + } + else if (frac == 1) { + FracGa = FracGa + TPCHit.Integral(); + } + else if (frac == 2) { + FracNe = FracNe + TPCHit.Integral(); + } + else { + FracRest = FracRest + TPCHit.Integral(); + } } } } + + long unsigned int GeneratorType = ProducerUtils::WhichGeneratorType(GeneratorParticles, MainTrID); + VecGenPur[int(GeneratorType)] = VecGenPur[int(GeneratorType)] + TPCHit.Integral(); + if (SignalTrackIDs.find(MainTrID) != SignalTrackIDs.end()) { + hitCharge = TPCHit.Integral(); + Pur = Pur + hitCharge; + } } - long unsigned int GeneratorType = ProducerUtils::WhichGeneratorType(GeneratorParticles, MainTrID); - VecGenPur[int(GeneratorType)] = VecGenPur[int(GeneratorType)] + TPCHit.Integral(); - if (SignalTrackIDs.find(MainTrID) != SignalTrackIDs.end()) { - hitCharge = TPCHit.Integral(); - Pur = Pur + hitCharge; + float MainGenPurity = 0; + for (size_t genpur = 0; genpur < VecGenPur.size(); genpur++) + { + VecGenPur[genpur] = VecGenPur[genpur] / ncharge; + if (VecGenPur[genpur] > MainGenPurity) { + MainGenerator = genpur; + MainGenPurity = VecGenPur[genpur]; + } } - } - float MainGenPurity = 0; - for (size_t genpur = 0; genpur < VecGenPur.size(); genpur++) - { - VecGenPur[genpur] = VecGenPur[genpur] / ncharge; - if (VecGenPur[genpur] > MainGenPurity) { - MainGenerator = genpur; - MainGenPurity = VecGenPur[genpur]; + dzdy = thisdzdy[0]; + thisdzdy.clear(); + FracE /= ncharge; + FracGa /= ncharge; + FracNe /= ncharge; + FracRest /= ncharge; + clustX /= ncharge; + clustY /= ncharge; + clustZ /= ncharge; + clustT /= ncharge; + Pur /= ncharge; + + ClIdxMap[ClustersIdx[idx][i]] = {idx, i}; // Map the cluster index to the plane and cluster number + ClNHits[idx].push_back(nhit); + ClCharge[idx].push_back(ncharge); + ClMaxCharge[idx].push_back(maxHit); + ClT[idx].push_back(clustT); + ClTPC[idx].push_back(clustTPC); + ClX[idx].push_back(clustX); + ClY[idx].push_back(clustY); + ClZ[idx].push_back(clustZ); + ClFracE[idx].push_back(FracE); + ClFracGa[idx].push_back(FracGa); + ClFracNe[idx].push_back(FracNe); + ClFracRest[idx].push_back(FracRest); + ClPur[idx].push_back(Pur); + ClGen[idx].push_back(MainGenerator); + ClGenPur[idx].push_back(MainGenPurity); + Cldzdy[idx].push_back(dzdy); + ClMainID[idx].push_back(MainTrID); + ClVecGenPur[idx].push_back(VecGenPur); + } + } // Finished first cluster processing + + //-------------------------------------------------------------------- Cluster Matching -------------------------------------------------------------------------// + std::vector MVecGen = {}; + std::vector> MVecGenFrac = {}; + std::vector MVecFracE = {}, MVecFracGa = {}, MVecFracNe = {}, MVecFracRest = {}, MVecGenPur = {}; + std::vector> MatchedClustersIdx = {{}, {}, {}}; + std::vector> MVecMainID = {{}, {}, {}}, MVecNHits = {{}, {}, {}}, MVecTPC = {{}, {}, {}}, MVecChannel = {{}, {}, {}}; + std::vector> MVecPur = {{}, {}, {}}, MVecMaxCharge = {{}, {}, {}}, MVecCharge = {{}, {}, {}}, MVecTime = {{}, {}, {}}, MVecRecoX = {{}, {}, {}}, MVecRecoY = {{}, {}, {}}, MVecRecoZ = {{}, {}, {}}; + std::vector> MVecDirDir = {{}, {}, {}}, MatchedClCompleteness = {{}, {}, {}}, MVecdT = {{}, {}, {}}; + std::vector SolarClusters; + std::vector> SolarClustersPtr; + + std::vector>> MatchedClusters = {{}, {}, {}}; + + // If present, grab the SolarClusters from the event + if (fGenerateSolarCluster == false) { + art::Handle> SolarClusterHandle; + evt.getByLabel(fSolarClusterLabel, SolarClusterHandle); + if (SolarClusterHandle.isValid()) { + for (size_t i = 0; i < SolarClusterHandle->size(); i++) { + SolarClustersPtr.push_back(art::Ptr(SolarClusterHandle, i)); } } + // Requires further implementation to match the previous "planecluster" + // ... + } + else { + std::string SolarClusterInfo = "SolarClusterInfo: "; + SolarClusterInfo = SolarClusterInfo + "(" + ProducerUtils::str(Clusters0.size()) + "," + ProducerUtils::str(Clusters1.size()) + "," + ProducerUtils::str(Clusters2.size()) + ")"; + lowe->MatchClusters(SignalTrackIDs, MatchedClustersIdx, MatchedClusters, ClustersIdx, AllPlaneClusters, MVecMainID, MVecNHits, MVecTPC, MVecChannel, MVecTime, MVecRecoY, MVecRecoZ, MVecDirDir, MVecCharge, MVecPur, MatchedClCompleteness, clockData, true); - dzdy = thisdzdy[0]; - thisdzdy.clear(); - FracE /= ncharge; - FracGa /= ncharge; - FracNe /= ncharge; - FracRest /= ncharge; - clustX /= ncharge; - clustY /= ncharge; - clustZ /= ncharge; - clustT /= ncharge; - Pur /= ncharge; - - ClIdxMap[ClustersIdx[idx][i]] = {idx, i}; // Map the cluster index to the plane and cluster number - ClNHits[idx].push_back(nhit); - ClCharge[idx].push_back(ncharge); - ClMaxCharge[idx].push_back(maxHit); - ClT[idx].push_back(clustT); - ClTPC[idx].push_back(clustTPC); - ClX[idx].push_back(clustX); - ClY[idx].push_back(clustY); - ClZ[idx].push_back(clustZ); - ClFracE[idx].push_back(FracE); - ClFracGa[idx].push_back(FracGa); - ClFracNe[idx].push_back(FracNe); - ClFracRest[idx].push_back(FracRest); - ClPur[idx].push_back(Pur); - ClGen[idx].push_back(MainGenerator); - ClGenPur[idx].push_back(MainGenPurity); - Cldzdy[idx].push_back(dzdy); - ClMainID[idx].push_back(MainTrID); - ClVecGenPur[idx].push_back(VecGenPur); - } - } // Finished first cluster processing - - //-------------------------------------------------------------------- Cluster Matching -------------------------------------------------------------------------// - std::vector MVecGen = {}; - std::vector> MVecGenFrac = {}; - std::vector MVecFracE = {}, MVecFracGa = {}, MVecFracNe = {}, MVecFracRest = {}, MVecGenPur = {}; - std::vector> MatchedClustersIdx = {{}, {}, {}}; - std::vector> MVecMainID = {{}, {}, {}}, MVecNHits = {{}, {}, {}}, MVecTPC = {{}, {}, {}}, MVecChannel = {{}, {}, {}}; - std::vector> MVecPur = {{}, {}, {}}, MVecMaxCharge = {{}, {}, {}}, MVecCharge = {{}, {}, {}}, MVecTime = {{}, {}, {}}, MVecRecoX = {{}, {}, {}}, MVecRecoY = {{}, {}, {}}, MVecRecoZ = {{}, {}, {}}; - std::vector> MVecDirDir = {{}, {}, {}}, MatchedClCompleteness = {{}, {}, {}}, MVecdT = {{}, {}, {}}; - std::vector SolarClusters; - std::vector> SolarClustersPtr; - - std::vector>> MatchedClusters = {{}, {}, {}}; - - // If present, grab the SolarClusters from the event - if (fGenerateSolarCluster == false) { - art::Handle> SolarClusterHandle; - evt.getByLabel(fSolarClusterLabel, SolarClusterHandle); - if (SolarClusterHandle.isValid()) { - for (size_t i = 0; i < SolarClusterHandle->size(); i++) { - SolarClustersPtr.push_back(art::Ptr(SolarClusterHandle, i)); + SolarClusterInfo = SolarClusterInfo + "\nFound " + ProducerUtils::str(int(MatchedClustersIdx[2].size())) + " MatchedClusters (from col. plane loop)!"; + for (int ThisClIdx = 0; ThisClIdx < int(MatchedClustersIdx[2].size()); ThisClIdx++) + { + for (int plane = 0; plane < 2; plane++) + { + if (MVecTime[plane][ThisClIdx] > -1e6) { // There is a matched cluster in this plane + int RefClIdx = ClIdxMap[MatchedClustersIdx[plane][ThisClIdx]][1]; // Get the cluster index in the plane + MVecMainID[plane].push_back(ClMainID[plane][RefClIdx]); + MVecdT[plane].push_back(abs(MVecTime[2][ThisClIdx] - MVecTime[plane][ThisClIdx])); + MVecMaxCharge[plane].push_back(ClMaxCharge[plane][RefClIdx]); + SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane " + ProducerUtils::str(plane) + " with time " + ProducerUtils::str(MVecTime[plane][ThisClIdx]) + " and charge " + ProducerUtils::str(MVecCharge[plane][ThisClIdx]) + " with TPC " + ProducerUtils::str(MVecTPC[plane][ThisClIdx]); + } + + else { // No matched cluster in this plane, fill with -1 + MVecMainID[plane].push_back(-1); + MVecdT[plane].push_back(-1e6); + MVecMaxCharge[plane].push_back(-1e6); + SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane " + ProducerUtils::str(plane) + " with time -1e6 and charge -1e6 with TPC -1"; + } + } + int RefClIdx = ClIdxMap[MatchedClustersIdx[2][ThisClIdx]][1]; // Get the plane index of the matched cluster + MVecMainID[2].push_back(ClMainID[2][RefClIdx]); + MVecRecoX[2].push_back(ClT[2][RefClIdx] *driftLength/driftTime); // Convert to microseconds and then to cm + MVecMaxCharge[2].push_back(ClMaxCharge[2][RefClIdx]); + MVecGenPur.push_back(ClGenPur[2][RefClIdx]); + MVecGen.push_back(ClGen[2][RefClIdx]); + MVecFracE.push_back(ClFracE[2][RefClIdx]); + MVecFracGa.push_back(ClFracGa[2][RefClIdx]); + MVecFracNe.push_back(ClFracNe[2][RefClIdx]); + MVecFracRest.push_back(ClFracRest[2][RefClIdx]); + MVecGenFrac.push_back(ClVecGenPur[2][RefClIdx]); + SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane 2 with time " + ProducerUtils::str(MVecTime[2][ThisClIdx]) + " and charge " + ProducerUtils::str(MVecCharge[2][ThisClIdx]) + " with TPC " + ProducerUtils::str(MVecTPC[2][ThisClIdx]); } - } - // Requires further implementation to match the previous "planecluster" - // ... - } - else { - std::string SolarClusterInfo = "SolarClusterInfo: "; - SolarClusterInfo = SolarClusterInfo + "(" + ProducerUtils::str(Clusters0.size()) + "," + ProducerUtils::str(Clusters1.size()) + "," + ProducerUtils::str(Clusters2.size()) + ")"; - lowe->MatchClusters(SignalTrackIDs, MatchedClustersIdx, MatchedClusters, ClustersIdx, AllPlaneClusters, MVecMainID, MVecNHits, MVecTPC, MVecChannel, MVecTime, MVecRecoY, MVecRecoZ, MVecDirDir, MVecCharge, MVecPur, MatchedClCompleteness, clockData, true); + producer->PrintInColor(SolarClusterInfo, ProducerUtils::GetColor("yellow"), "Debug"); - 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++) + for (int i = 0; i < int(MVecNHits[2].size()); i++) { - if (MVecTime[plane][ThisClIdx] > -1e6) { // There is a matched cluster in this plane - int RefClIdx = ClIdxMap[MatchedClustersIdx[plane][ThisClIdx]][1]; // Get the cluster index in the plane - MVecMainID[plane].push_back(ClMainID[plane][RefClIdx]); - MVecdT[plane].push_back(abs(MVecTime[2][ThisClIdx] - MVecTime[plane][ThisClIdx])); - MVecMaxCharge[plane].push_back(ClMaxCharge[plane][RefClIdx]); - SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane " + ProducerUtils::str(plane) + " with time " + ProducerUtils::str(MVecTime[plane][ThisClIdx]) + " and charge " + ProducerUtils::str(MVecCharge[plane][ThisClIdx]) + " with TPC " + ProducerUtils::str(MVecTPC[plane][ThisClIdx]); + if (fClusterPreselectionSignal && MVecPur[2][i] == 0) + { + continue; } - - else { // No matched cluster in this plane, fill with -1 - MVecMainID[plane].push_back(-1); - MVecdT[plane].push_back(-1e6); - MVecMaxCharge[plane].push_back(-1e6); - SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane " + ProducerUtils::str(plane) + " with time -1e6 and charge -1e6 with TPC -1"; + std::vector clustPos = {MVecRecoX[2][i], MVecRecoY[2][i], MVecRecoZ[2][i]}; + int clustMainID = MVecMainID[2][i]; + int clustNHits = MVecNHits[2][i]; + int clustTPC = MVecTPC[2][i]; + int clustChannel = MVecChannel[2][i]; + float clustCharge = MVecCharge[2][i]; + float clustTime = MVecTime[2][i]; + float clustPurity = MVecPur[2][i]; + float clustCompleteness = MatchedClCompleteness[2][i]; + std::vector clustVector = {}; // Vector of recob::Cluster + // Add clusters according to the indices in MatchedClustersIdx + for (int plane = 0; plane < 3; plane++) + { + int clustIdx = ClIdxMap[MatchedClustersIdx[plane][i]][1]; + if (clustIdx >= 0 && clustIdx < int(AllPlaneClusters[plane].size())) + { + // Cluster(float start_wire,float sigma_start_wire,float start_tick,float sigma_start_tick,float start_charge,float start_angle,float start_opening,float end_wire,float sigma_end_wire,float end_tick,float sigma_end_tick,float end_charge,float end_angle,float end_opening,float integral,float integral_stddev,float summedADC,float summedADC_stddev,unsigned int n_hits,float multiple_hit_density,float width,ID_t ID,geo::View_t view,geo::PlaneID const& plane,SentryArgument_t sentry = Sentry); + recob::Cluster thisCluster( + ClY[plane][clustIdx], 0, ClT[plane][clustIdx], 0, + ClCharge[plane][clustIdx], 0, 0, + ClY[plane][clustIdx], 0, ClT[plane][clustIdx], 0, + ClCharge[plane][clustIdx], 0, 0, + ClCharge[plane][clustIdx], 0, + ClCharge[plane][clustIdx], 0, + int(AllPlaneClusters[plane][clustIdx].size()), + 0, + 0, + clustIdx, + geo::View_t(plane), + geo::PlaneID(0, 0, plane), + {} // Assuming default for SentryArgument_t + ); + clustVector.push_back(thisCluster); + } } - } - int RefClIdx = ClIdxMap[MatchedClustersIdx[2][ThisClIdx]][1]; // Get the plane index of the matched cluster - MVecMainID[2].push_back(ClMainID[2][RefClIdx]); - MVecRecoX[2].push_back(ClT[2][RefClIdx] *driftLength/driftTime); // Convert to microseconds and then to cm - MVecMaxCharge[2].push_back(ClMaxCharge[2][RefClIdx]); - MVecGenPur.push_back(ClGenPur[2][RefClIdx]); - MVecGen.push_back(ClGen[2][RefClIdx]); - MVecFracE.push_back(ClFracE[2][RefClIdx]); - MVecFracGa.push_back(ClFracGa[2][RefClIdx]); - MVecFracNe.push_back(ClFracNe[2][RefClIdx]); - MVecFracRest.push_back(ClFracRest[2][RefClIdx]); - MVecGenFrac.push_back(ClVecGenPur[2][RefClIdx]); - SolarClusterInfo = SolarClusterInfo + "\nMatched Cluster in plane 2 with time " + ProducerUtils::str(MVecTime[2][ThisClIdx]) + " and charge " + ProducerUtils::str(MVecCharge[2][ThisClIdx]) + " with TPC " + ProducerUtils::str(MVecTPC[2][ThisClIdx]); + + solar::LowECluster ThisSolarCluster(clustPos, clustMainID, clustNHits, clustTPC, clustChannel, clustCharge, clustTime, clustPurity, clustCompleteness, clustVector); + SolarClusters.push_back(ThisSolarCluster); + } } - producer->PrintInColor(SolarClusterInfo, ProducerUtils::GetColor("yellow"), "Debug"); + //-------------------------------------------------------------------- Cluster Tree Export -------------------------------------------------------------------------// + // Need to implement the primary cluster finding based on external algorithm that uses the solar::LowECluster + // std::vector EventCandidateFound = {}; + // std::vector>> EventCandidateVector; + // std::vector> EventCandidateIdx; + // lowe->FindPrimaryClusters(SolarClustersPtr, EventCandidateFound, EventCandidateVector, EventCandidateIdx, clockData, evt); + + // For now, loop over matched clusters and export to tree if all conditions are satisfied + std::string sClustersReco = "\n# ClusterReco: Looping over " + ProducerUtils::str(int(MVecNHits[2].size())) + " matched clusters"; + producer->PrintInColor(sClustersReco, ProducerUtils::GetColor("green")); for (int i = 0; i < int(MVecNHits[2].size()); i++) { if (fClusterPreselectionSignal && MVecPur[2][i] == 0) { continue; } - std::vector clustPos = {MVecRecoX[2][i], MVecRecoY[2][i], MVecRecoZ[2][i]}; - int clustMainID = MVecMainID[2][i]; - int clustNHits = MVecNHits[2][i]; - int clustTPC = MVecTPC[2][i]; - int clustChannel = MVecChannel[2][i]; - float clustCharge = MVecCharge[2][i]; - float clustTime = MVecTime[2][i]; - float clustPurity = MVecPur[2][i]; - float clustCompleteness = MatchedClCompleteness[2][i]; - std::vector clustVector = {}; // Vector of recob::Cluster - // Add clusters according to the indices in MatchedClustersIdx - for (int plane = 0; plane < 3; plane++) + bool TrackMatch = false; + bool AdjClusterMatch = false; + std::string sFlashReco = ""; + std::string sVertexReco = ""; + std::string sClusterReco = ""; + std::string sResultColor = "white"; + std::string sAdjClusters = ""; + float OpFlashResidual = 0; + float MatchedOpFlashPE = -1e6; + // float MatchedOpFlashResidual = 1e6; + float MatchedOpFlashX = -1e6; + + if (MVecNHits[2][i] > fClusterPreselectionNHits) { - int clustIdx = ClIdxMap[MatchedClustersIdx[plane][i]][1]; - if (clustIdx >= 0 && clustIdx < int(AllPlaneClusters[plane].size())) + MPrimary = true; + MAdjClNum = 0; + MSignalAdjClNum = 0; + MAdjClTime = {}; + MAdjClCharge = {}; + MAdjClInd0Charge = {}; + MAdjClInd1Charge = {}; + MAdjClMaxCharge = {}; + MAdjClInd0MaxCharge = {}; + MAdjClInd1MaxCharge = {}; + MAdjClNHits = {}; + MAdjClInd0NHits = {}; + MAdjClInd1NHits = {}; + MAdjClRecoY = {}; + MAdjClRecoZ = {}; + MAdjClR = {}; + MAdjClPur = {}; + MAdjClGen = {}; + MAdjClGenPur = {}; + MAdjClMainID = {}; + MAdjClMainPDG = {}; + MAdjClMainE = {}; + MAdjClMainP = {}; + MAdjClMainK = {}; + MAdjClMainX = {}; + MAdjClMainY = {}; + MAdjClMainZ = {}; + MAdjClEndX = {}; + MAdjClEndY = {}; + MAdjClEndZ = {}; + MAdjFlashR = {}; + MAdjFlashPE = {}; + MAdjFlashPur = {}; + MAdjFlashSTD = {}; + MAdjFlashTime = {}; + MAdjFlashNHits = {}; + MAdjFlashPlane = {}; + MAdjFlashFast = {}; + MAdjFlashMaxPE = {}; + MAdjFlashRecoX = {}; + MAdjFlashRecoY = {}; + MAdjFlashRecoZ = {}; + MAdjFlashResidual = {}; + MTrackStart = {-1e6, -1e6, -1e6}; + MTrackEnd = {-1e6, -1e6, -1e6}; + + for (int j = 0; j < int(MVecNHits[2].size()); j++) { - // Cluster(float start_wire,float sigma_start_wire,float start_tick,float sigma_start_tick,float start_charge,float start_angle,float start_opening,float end_wire,float sigma_end_wire,float end_tick,float sigma_end_tick,float end_charge,float end_angle,float end_opening,float integral,float integral_stddev,float summedADC,float summedADC_stddev,unsigned int n_hits,float multiple_hit_density,float width,ID_t ID,geo::View_t view,geo::PlaneID const& plane,SentryArgument_t sentry = Sentry); - recob::Cluster thisCluster( - ClY[plane][clustIdx], 0, ClT[plane][clustIdx], 0, - ClCharge[plane][clustIdx], 0, 0, - ClY[plane][clustIdx], 0, ClT[plane][clustIdx], 0, - ClCharge[plane][clustIdx], 0, 0, - ClCharge[plane][clustIdx], 0, - ClCharge[plane][clustIdx], 0, - int(AllPlaneClusters[plane][clustIdx].size()), - 0, - 0, - clustIdx, - geo::View_t(plane), - geo::PlaneID(0, 0, plane), - {} // Assuming default for SentryArgument_t - ); - clustVector.push_back(thisCluster); - } - } - - solar::LowECluster ThisSolarCluster(clustPos, clustMainID, clustNHits, clustTPC, clustChannel, clustCharge, clustTime, clustPurity, clustCompleteness, clustVector); - SolarClusters.push_back(ThisSolarCluster); - } - } - - //-------------------------------------------------------------------- Cluster Tree Export -------------------------------------------------------------------------// - // Need to implement the primary cluster finding based on external algorithm that uses the solar::LowECluster - // std::vector EventCandidateFound = {}; - // std::vector>> EventCandidateVector; - // std::vector> EventCandidateIdx; - // lowe->FindPrimaryClusters(SolarClustersPtr, EventCandidateFound, EventCandidateVector, EventCandidateIdx, clockData, evt); - - // For now, loop over matched clusters and export to tree if all conditions are satisfied - std::string sClustersReco = "\n# ClusterReco: Looping over " + ProducerUtils::str(int(MVecNHits[2].size())) + " matched clusters"; - producer->PrintInColor(sClustersReco, ProducerUtils::GetColor("green")); - for (int i = 0; i < int(MVecNHits[2].size()); i++) - { - if (fClusterPreselectionSignal && MVecPur[2][i] == 0) - { - continue; - } - bool TrackMatch = false; - bool AdjClusterMatch = false; - std::string sFlashReco = ""; - std::string sVertexReco = ""; - std::string sClusterReco = ""; - std::string sResultColor = "white"; - std::string sAdjClusters = ""; - float OpFlashResidual = 0; - float MatchedOpFlashPE = -1e6; - float MatchedOpFlashResidual = 1e6; - float MatchedOpFlashX = -1e6; - - if (MVecNHits[2][i] > fClusterPreselectionNHits) - { - MPrimary = true; - MAdjClNum = 0; - MSignalAdjClNum = 0; - MAdjClTime = {}; - MAdjClCharge = {}; - MAdjClInd0Charge = {}; - MAdjClInd1Charge = {}; - MAdjClMaxCharge = {}; - MAdjClInd0MaxCharge = {}; - MAdjClInd1MaxCharge = {}; - MAdjClNHits = {}; - MAdjClInd0NHits = {}; - MAdjClInd1NHits = {}; - MAdjClRecoY = {}; - MAdjClRecoZ = {}; - MAdjClR = {}; - MAdjClPur = {}; - MAdjClGen = {}; - MAdjClGenPur = {}; - MAdjClMainID = {}; - MAdjClMainPDG = {}; - MAdjClMainE = {}; - MAdjClMainP = {}; - MAdjClMainK = {}; - MAdjClMainX = {}; - MAdjClMainY = {}; - MAdjClMainZ = {}; - MAdjClEndX = {}; - MAdjClEndY = {}; - MAdjClEndZ = {}; - MAdjFlashR = {}; - MAdjFlashPE = {}; - MAdjFlashPur = {}; - MAdjFlashSTD = {}; - MAdjFlashTime = {}; - MAdjFlashNHits = {}; - MAdjFlashPlane = {}; - MAdjFlashFast = {}; - MAdjFlashMaxPE = {}; - MAdjFlashRecoX = {}; - MAdjFlashRecoY = {}; - MAdjFlashRecoZ = {}; - MAdjFlashResidual = {}; - MTrackStart = {-1e6, -1e6, -1e6}; - MTrackEnd = {-1e6, -1e6, -1e6}; - - for (int j = 0; j < int(MVecNHits[2].size()); j++) - { - if (j == i) { continue;} // Do not compare the cluster with itself + if (j == i) { continue;} // Do not compare the cluster with itself - double ClusterDistance = 0; - producer->ComputeDistance3D(ClusterDistance, MVecTime[2][i], MVecRecoY[2][i], MVecRecoZ[2][i], MVecTime[2][j], MVecRecoY[2][j], MVecRecoZ[2][j], TPCIDdriftLength[MVecTPC[2][i]], TPCIDdriftTime[MVecTPC[2][i]]); - if (MVecCharge[2][j] < fMinClusterCharge) { continue; } // Skip clusters that are too small - if (MVecTPC[2][j]%2 != MVecTPC[2][i]%2) { continue; } // Skip clusters that are in different sides of the detector - if (ClusterDistance > fAdjClusterRad) { continue; } // Skip clusters that are too far + double ClusterDistance = 0; + producer->ComputeDistance3D(ClusterDistance, MVecTime[2][i], MVecRecoY[2][i], MVecRecoZ[2][i], MVecTime[2][j], MVecRecoY[2][j], MVecRecoZ[2][j], TPCIDdriftLength[MVecTPC[2][i]], TPCIDdriftTime[MVecTPC[2][i]]); + if (MVecCharge[2][j] < fMinClusterCharge) { continue; } // Skip clusters that are too small + if (MVecTPC[2][j]%2 != MVecTPC[2][i]%2) { continue; } // Skip clusters that are in different sides of the detector + if (ClusterDistance > fAdjClusterRad) { continue; } // Skip clusters that are too far - sAdjClusters += " - Cluster " + ProducerUtils::str(j) + " at distance " + ProducerUtils::str(ClusterDistance) + " with time " + ProducerUtils::str(MVecTime[2][j]) + " and charge " + ProducerUtils::str(MVecCharge[2][j]) + " in TPC " + ProducerUtils::str(MVecTPC[2][j]); - sAdjClusters += " and hits " + ProducerUtils::str(MVecNHits[2][j]) + "\n"; + sAdjClusters += " - Cluster " + ProducerUtils::str(j) + " at distance " + ProducerUtils::str(ClusterDistance) + " with time " + ProducerUtils::str(MVecTime[2][j]) + " and charge " + ProducerUtils::str(MVecCharge[2][j]) + " in TPC " + ProducerUtils::str(MVecTPC[2][j]); + sAdjClusters += " and hits " + ProducerUtils::str(MVecNHits[2][j]) + "\n"; - if (MVecCharge[2][j] > MVecCharge[2][i]) { MPrimary = false; } - if (MVecGen[i] == MVecGen[j]) { MSignalAdjClNum += 1; } - MAdjClNum += 1; - - // If the cluster is matched, add the information to the vectors - AdjClusterMatch = true; - MAdjClTime.push_back(MVecTime[2][j]); - MAdjClInd0Charge.push_back(MVecCharge[0][j]); - MAdjClInd1Charge.push_back(MVecCharge[1][j]); - MAdjClCharge.push_back(MVecCharge[2][j]); - MAdjClMaxCharge.push_back(MVecMaxCharge[2][j]); - MAdjClInd0MaxCharge.push_back(MVecMaxCharge[0][j]); - MAdjClInd1MaxCharge.push_back(MVecMaxCharge[1][j]); - MAdjClNHits.push_back(MVecNHits[2][j]); - MAdjClInd0NHits.push_back(MVecNHits[0][j]); - MAdjClInd1NHits.push_back(MVecNHits[1][j]); - MAdjClRecoY.push_back(MVecRecoY[2][j]); - MAdjClRecoZ.push_back(MVecRecoZ[2][j]); - MAdjClR.push_back(sqrt(pow(MVecRecoY[2][i] - MVecRecoY[2][j], 2) + pow(MVecRecoZ[2][i] - MVecRecoZ[2][j], 2))); - MAdjClPur.push_back(MVecPur[2][j]); - MAdjClGen.push_back(MVecGen[j]); - MAdjClGenPur.push_back(MVecGenPur[j]); - MAdjClMainID.push_back(MVecMainID[2][j]); - - // If mother exists add the mother information - const simb::MCParticle *MAdjClTruth; - int TerminalOutput = ProducerUtils::supress_stdout(); - MAdjClTruth = pi_serv->TrackIdToParticle_P(MVecMainID[2][j]); - ProducerUtils::resume_stdout(TerminalOutput); - - if (MAdjClTruth == 0) { - MAdjClMainPDG.push_back(0); - MAdjClMainE.push_back(-1e6); - MAdjClMainP.push_back(-1e6); - MAdjClMainK.push_back(-1e6); - MAdjClMainX.push_back(-1e6); - MAdjClMainY.push_back(-1e6); - MAdjClMainZ.push_back(-1e6); - MAdjClEndX.push_back(-1e6); - MAdjClEndY.push_back(-1e6); - MAdjClEndZ.push_back(-1e6); - } - else { - MAdjClMainPDG.push_back(MAdjClTruth->PdgCode()); - MAdjClMainE.push_back(1e3*MAdjClTruth->E()); - MAdjClMainP.push_back(1e3*MAdjClTruth->P()); - MAdjClMainK.push_back(1e3*MAdjClTruth->E() - 1e3*MAdjClTruth->Mass()); - MAdjClMainX.push_back(MAdjClTruth->Vx()); - MAdjClMainY.push_back(MAdjClTruth->Vy()); - MAdjClMainZ.push_back(MAdjClTruth->Vz()); - MAdjClEndX.push_back(MAdjClTruth->EndX()); - MAdjClEndY.push_back(MAdjClTruth->EndY()); - MAdjClEndZ.push_back(MAdjClTruth->EndZ()); + if (MVecCharge[2][j] > MVecCharge[2][i]) { MPrimary = false; } + if (MVecGen[i] == MVecGen[j]) { MSignalAdjClNum += 1; } + MAdjClNum += 1; + + // If the cluster is matched, add the information to the vectors + AdjClusterMatch = true; + MAdjClTime.push_back(MVecTime[2][j]); + MAdjClInd0Charge.push_back(MVecCharge[0][j]); + MAdjClInd1Charge.push_back(MVecCharge[1][j]); + MAdjClCharge.push_back(MVecCharge[2][j]); + MAdjClMaxCharge.push_back(MVecMaxCharge[2][j]); + MAdjClInd0MaxCharge.push_back(MVecMaxCharge[0][j]); + MAdjClInd1MaxCharge.push_back(MVecMaxCharge[1][j]); + MAdjClNHits.push_back(MVecNHits[2][j]); + MAdjClInd0NHits.push_back(MVecNHits[0][j]); + MAdjClInd1NHits.push_back(MVecNHits[1][j]); + MAdjClRecoY.push_back(MVecRecoY[2][j]); + MAdjClRecoZ.push_back(MVecRecoZ[2][j]); + MAdjClR.push_back(sqrt(pow(MVecRecoY[2][i] - MVecRecoY[2][j], 2) + pow(MVecRecoZ[2][i] - MVecRecoZ[2][j], 2))); + MAdjClPur.push_back(MVecPur[2][j]); + MAdjClGen.push_back(MVecGen[j]); + MAdjClGenPur.push_back(MVecGenPur[j]); + MAdjClMainID.push_back(MVecMainID[2][j]); + + // If mother exists add the mother information + const simb::MCParticle *MAdjClTruth; + int TerminalOutput = ProducerUtils::supress_stdout(); + MAdjClTruth = pi_serv->TrackIdToParticle_P(MVecMainID[2][j]); + ProducerUtils::resume_stdout(TerminalOutput); + + if (MAdjClTruth == 0) { + MAdjClMainPDG.push_back(0); + MAdjClMainE.push_back(-1e6); + MAdjClMainP.push_back(-1e6); + MAdjClMainK.push_back(-1e6); + MAdjClMainX.push_back(-1e6); + MAdjClMainY.push_back(-1e6); + MAdjClMainZ.push_back(-1e6); + MAdjClEndX.push_back(-1e6); + MAdjClEndY.push_back(-1e6); + MAdjClEndZ.push_back(-1e6); + } + else { + MAdjClMainPDG.push_back(MAdjClTruth->PdgCode()); + MAdjClMainE.push_back(1e3*MAdjClTruth->E()); + MAdjClMainP.push_back(1e3*MAdjClTruth->P()); + MAdjClMainK.push_back(1e3*MAdjClTruth->E() - 1e3*MAdjClTruth->Mass()); + MAdjClMainX.push_back(MAdjClTruth->Vx()); + MAdjClMainY.push_back(MAdjClTruth->Vy()); + MAdjClMainZ.push_back(MAdjClTruth->Vz()); + MAdjClEndX.push_back(MAdjClTruth->EndX()); + MAdjClEndY.push_back(MAdjClTruth->EndY()); + MAdjClEndZ.push_back(MAdjClTruth->EndZ()); + } } - } - sResultColor = "yellow"; - if (MVecPur[2][i] > 0) { - sResultColor = "green"; - } + sResultColor = "yellow"; + if (MVecPur[2][i] > 0) { + sResultColor = "green"; + } - if (fClusterPreselectionPrimary && !MPrimary) { continue; } + if (fClusterPreselectionPrimary && !MPrimary) { continue; } - if (MPrimary) { - sClusterReco += "*** Matched preselection cluster: " + ProducerUtils::str(i) + "\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]; - } - else { - sClusterReco += " - Gen ?? -> Unknown"; - } + if (MPrimary) { + 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]; + } + else { + sClusterReco += " - Gen ?? -> Unknown"; + } - sClusterReco += " TPC " + ProducerUtils::str(MVecTPC[2][i]) + "\n"; - sClusterReco += " - Purity " + ProducerUtils::str(MVecGenPur[i]) + " Hits " + ProducerUtils::str(MVecNHits[2][i]) + "\n"; - sClusterReco += " - Charge " + ProducerUtils::str(MVecCharge[2][i]) + " ( MaxHit " + ProducerUtils::str(MVecMaxCharge[2][i]) + " )\n"; - sClusterReco += " - #AdjCl " + ProducerUtils::str(MAdjClNum) + " ( " + ProducerUtils::str(MSignalAdjClNum) + " signal ):\n"; + sClusterReco += " TPC " + ProducerUtils::str(MVecTPC[2][i]) + "\n"; + sClusterReco += " - Purity " + ProducerUtils::str(MVecGenPur[i]) + " Hits " + ProducerUtils::str(MVecNHits[2][i]) + "\n"; + sClusterReco += " - Charge " + ProducerUtils::str(MVecCharge[2][i]) + " ( MaxHit " + ProducerUtils::str(MVecMaxCharge[2][i]) + " )\n"; + sClusterReco += " - #AdjCl " + ProducerUtils::str(MAdjClNum) + " ( " + ProducerUtils::str(MSignalAdjClNum) + " signal ):\n"; + + if (AdjClusterMatch) { sClusterReco += sAdjClusters; } + sClusterReco += " - RecoCol Time,Y,Z ( " + ProducerUtils::str(MVecTime[2][i]) + ", " + ProducerUtils::str(MVecRecoY[2][i]) + ", " + ProducerUtils::str(MVecRecoZ[2][i]) + " )\n"; + sClusterReco += " - RecoInd0 Time,Y,Z ( " + ProducerUtils::str(MVecTime[0][i]) + ", " + ProducerUtils::str(MVecRecoY[0][i]) + ", " + ProducerUtils::str(MVecRecoZ[0][i]) + " )\n"; + sClusterReco += " - RecoInd1 Time,Y,Z ( " + ProducerUtils::str(MVecTime[1][i]) + ", " + ProducerUtils::str(MVecRecoY[1][i]) + ", " + ProducerUtils::str(MVecRecoZ[1][i]) + " )\n"; + + if (fSaveTrackInfo) { + TVector3 ThisClVertex = {0, MVecRecoY[2][i], MVecRecoZ[2][i]}; + float MaxVertexDistance = 10; // if track is further away from ThisClVertex than + for (int i = 0; i < TrackNum; i++) + { // using index loop to get track idx + recob::Track trk = *TrackList[i]; + TVector3 trk_start(0, trk.Start().Y(), trk.Start().Z()); + TVector3 trk_end(0, trk.End().Y(), trk.End().Z()); + // throw away bad tracks + if ((trk_start - ThisClVertex).Mag() > MaxVertexDistance && (trk_end - ThisClVertex).Mag() > MaxVertexDistance) { continue; } + + MTrackNPoints = trk.NPoints(); + MTrackStart = {trk.Start().X(), trk.Start().Y(), trk.Start().Z()}; + MTrackEnd = {trk.End().X(), trk.End().Y(), trk.End().Z()}; + MTrackChi2 = trk.Chi2(); + + sClusterReco += "*** Matched pmtrack: \n"; + sClusterReco += " - Track has start ( " + ProducerUtils::str(trk.Start().X()) + ", " + ProducerUtils::str(trk.Start().Y()) + ", " + ProducerUtils::str(trk.Start().Z()) + " )\n"; + sClusterReco += " - Track has end ( " + ProducerUtils::str(trk.End().X()) + ", " + ProducerUtils::str(trk.End().Y()) + ", " + ProducerUtils::str(trk.End().Z()) + " )\n\n"; + TrackMatch = true; + }; // Loop over tracks + }; // if (fSaveTrackInfo) + }; // if (MPrimary) - if (AdjClusterMatch) { sClusterReco += sAdjClusters; } - sClusterReco += " - RecoCol Time,Y,Z ( " + ProducerUtils::str(MVecTime[2][i]) + ", " + ProducerUtils::str(MVecRecoY[2][i]) + ", " + ProducerUtils::str(MVecRecoZ[2][i]) + " )\n"; - sClusterReco += " - RecoInd0 Time,Y,Z ( " + ProducerUtils::str(MVecTime[0][i]) + ", " + ProducerUtils::str(MVecRecoY[0][i]) + ", " + ProducerUtils::str(MVecRecoZ[0][i]) + " )\n"; - sClusterReco += " - RecoInd1 Time,Y,Z ( " + ProducerUtils::str(MVecTime[1][i]) + ", " + ProducerUtils::str(MVecRecoY[1][i]) + ", " + ProducerUtils::str(MVecRecoZ[1][i]) + " )\n"; + if (fClusterPreselectionTrack && !TrackMatch) { continue; } - if (fSaveTrackInfo) { - TVector3 ThisClVertex = {0, MVecRecoY[2][i], MVecRecoZ[2][i]}; - float MaxVertexDistance = 10; // if track is further away from ThisClVertex than - for (int i = 0; i < TrackNum; i++) - { // using index loop to get track idx - recob::Track trk = *TrackList[i]; - TVector3 trk_start(0, trk.Start().Y(), trk.Start().Z()); - TVector3 trk_end(0, trk.End().Y(), trk.End().Z()); - // throw away bad tracks - if ((trk_start - ThisClVertex).Mag() > MaxVertexDistance && (trk_end - ThisClVertex).Mag() > MaxVertexDistance) { continue; } - - MTrackNPoints = trk.NPoints(); - MTrackStart = {trk.Start().X(), trk.Start().Y(), trk.Start().Z()}; - MTrackEnd = {trk.End().X(), trk.End().Y(), trk.End().Z()}; - MTrackChi2 = trk.Chi2(); - - sClusterReco += "*** Matched pmtrack: \n"; - sClusterReco += " - Track has start ( " + ProducerUtils::str(trk.Start().X()) + ", " + ProducerUtils::str(trk.Start().Y()) + ", " + ProducerUtils::str(trk.Start().Z()) + " )\n"; - sClusterReco += " - Track has end ( " + ProducerUtils::str(trk.End().X()) + ", " + ProducerUtils::str(trk.End().Y()) + ", " + ProducerUtils::str(trk.End().Z()) + " )\n\n"; - TrackMatch = true; - }; // Loop over tracks - }; // if (fSaveTrackInfo) - }; // if (MPrimary) - - if (fClusterPreselectionTrack && !TrackMatch) { continue; } - - std::string sFlashMatching = ""; - for (int j = 0; j < int(OpFlashPE.size()); j++) - { - // Skip flashes with time outside the cluster time window - float OpFlashR = -1e6; - double MAdjFlashX = 0; - if ((MVecTime[2][i] - OpFlashTime[j]) < 0 || (MVecTime[2][i] - OpFlashTime[j]) > TPCIDdriftTime[MVecTPC[2][i]]) { continue; } - producer->ComputeDistanceX(MAdjFlashX, MVecTime[2][i], OpFlashTime[j], TPCIDdriftLength[MVecTPC[2][i]], TPCIDdriftTime[MVecTPC[2][i]]); - if (fGeometry == "HD" && MVecTPC[2][i]%2 == 0) { - MAdjFlashX = -MAdjFlashX; - } - if (fGeometry == "VD") { - MAdjFlashX = TPCIDdriftLength[MVecTPC[2][i]] / 2 - MAdjFlashX; - } + std::string sFlashMatching = ""; + bool IsFirstFlash = true; + for (int j = 0; j < int(OpFlashPE.size()); j++) + { + // Skip flashes with time outside the cluster time window + float OpFlashR = -1e6; + double MAdjFlashX = 0; + if ((MVecTime[2][i] - OpFlashTime[j]) < 0 || (MVecTime[2][i] - OpFlashTime[j]) > TPCIDdriftTime[MVecTPC[2][i]]) { continue; } + producer->ComputeDistanceX(MAdjFlashX, MVecTime[2][i], OpFlashTime[j], TPCIDdriftLength[MVecTPC[2][i]], TPCIDdriftTime[MVecTPC[2][i]]); + if (fGeometry == "HD" && MVecTPC[2][i]%2 == 0) { + MAdjFlashX = -MAdjFlashX; + } + if (fGeometry == "VD") { + MAdjFlashX = TPCIDdriftLength[MVecTPC[2][i]] / 2 - MAdjFlashX; + } - // Make an eliptical cut on the flash position based on the clusters plane - if (fGeometry == "HD") { - OpFlashR = sqrt(pow(MVecRecoY[2][i] - OpFlashY[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); - if (pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { - sFlashMatching += "Skipping flash " + ProducerUtils::str(j) + " at (X,Y,Z) = (" + ProducerUtils::str(OpFlashX[j]) + "," + ProducerUtils::str(OpFlashY[j]) + "," + ProducerUtils::str(OpFlashZ[j]) + ") outside cut at (X,Y,Z) = (" + ProducerUtils::str(MAdjFlashX) + "," + ProducerUtils::str(MVecRecoY[2][i]) + "," + ProducerUtils::str(MVecRecoZ[2][i]) + ") with R = " + ProducerUtils::str(OpFlashR) + " cm\n"; - continue; + // Make an eliptical cut on the flash position based on the clusters plane + if (fGeometry == "HD") { + OpFlashR = sqrt(pow(MVecRecoY[2][i] - OpFlashY[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); + if (pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { + sFlashMatching += "Skipping flash " + ProducerUtils::str(j) + " at (X,Y,Z) = (" + ProducerUtils::str(OpFlashX[j]) + "," + ProducerUtils::str(OpFlashY[j]) + "," + ProducerUtils::str(OpFlashZ[j]) + ") outside cut at (X,Y,Z) = (" + ProducerUtils::str(MAdjFlashX) + "," + ProducerUtils::str(MVecRecoY[2][i]) + "," + ProducerUtils::str(MVecRecoZ[2][i]) + ") with R = " + ProducerUtils::str(OpFlashR) + " cm\n"; + continue; + } } - } - else if (fGeometry == "VD" && OpFlashPlane[j] == 0) { // Cathode flashes - if (pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } - OpFlashR = sqrt(pow(MVecRecoY[2][i] - OpFlashY[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); - } - else if (fGeometry == "VD" && (OpFlashPlane[j] == 1 || OpFlashPlane[j] == 2)) { // Membrane flashes - if (fAdjOpFlashMembraneProjection) { - if (MVecRecoY[2][i] * OpFlashY[j] < 0) { continue; } // Only consider clusters and flashes on the same side of the detector - if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } + else if (fGeometry == "VD" && OpFlashPlane[j] == 0) { // Cathode flashes + if (pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } + OpFlashR = sqrt(pow(MVecRecoY[2][i] - OpFlashY[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); } - else { - if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } + else if (fGeometry == "VD" && (OpFlashPlane[j] == 1 || OpFlashPlane[j] == 2)) { // Membrane flashes + if (fAdjOpFlashMembraneProjection) { + if (MVecRecoY[2][i] * OpFlashY[j] < 0) { continue; } // Only consider clusters and flashes on the same side of the detector + if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } + } + else { + if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } + } + OpFlashR = sqrt(pow(MAdjFlashX - OpFlashX[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); + } + else if (fGeometry == "VD" && (OpFlashPlane[j] == 3 || OpFlashPlane[j] == 4)) { // End-Cap flashes + if (fAdjOpFlashEndCapProjection){ + if (MVecRecoZ[2][i] < fidVolZ / 2 && OpFlashPlane[j] == 3) { continue; } // Only consider clusters and flashes on the same half of the volume + if (MVecRecoZ[2][i] > fidVolZ / 2 && OpFlashPlane[j] == 4) { continue; } // Only consider clusters and flashes on the same half of the volume + if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) > 1){ continue; } + } + else{ + if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1){ + continue; + } + } + OpFlashR = sqrt(pow(MAdjFlashX - OpFlashX[j], 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2)); } - OpFlashR = sqrt(pow(MAdjFlashX - OpFlashX[j], 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2)); - } - else if (fGeometry == "VD" && (OpFlashPlane[j] == 3 || OpFlashPlane[j] == 4)) { // End-Cap flashes - if (fAdjOpFlashEndCapProjection){ - if (MVecRecoZ[2][i] < fidVolZ / 2 && OpFlashPlane[j] == 3) { continue; } // Only consider clusters and flashes on the same half of the volume - if (MVecRecoZ[2][i] > fidVolZ / 2 && OpFlashPlane[j] == 4) { continue; } // Only consider clusters and flashes on the same half of the volume - if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) > 1){ continue; } + else if (fGeometry == "VD" && OpFlashPlane[j] == -1) { + sFlashMatching += "Skipping flash " + ProducerUtils::str(j) + " with unknown plane " + ProducerUtils::str(OpFlashPlane[j]) + "\n"; + continue; } - else{ - if (pow(MAdjFlashX - OpFlashX[j], 2) / pow(fAdjOpFlashX, 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2) / pow(fAdjOpFlashY, 2) + pow(MVecRecoZ[2][i] - OpFlashZ[j], 2) / pow(fAdjOpFlashZ, 2) > 1){ - continue; + + MAdjFlashR.push_back(OpFlashR); + MAdjFlashPE.push_back(OpFlashPE[j]); + MAdjFlashTime.push_back(OpFlashTime[j]); + MAdjFlashNHits.push_back(OpFlashNHits[j]); + MAdjFlashPlane.push_back(OpFlashPlane[j]); + MAdjFlashMaxPE.push_back(OpFlashMaxPE[j]); + MAdjFlashFast.push_back(OpFlashFast[j]); + MAdjFlashRecoX.push_back(OpFlashX[j]); + MAdjFlashRecoY.push_back(OpFlashY[j]); + MAdjFlashRecoZ.push_back(OpFlashZ[j]); + MAdjFlashSTD.push_back(OpFlashSTD[j]); + MAdjFlashPur.push_back(OpFlashPur[j]); + + // Compute the residual between the predicted cluster signal and the flash + adjophits->FlashMatchResidual( OpFlashResidual, OpHitVec[j], MAdjFlashX, double(MVecRecoY[2][i]), double(MVecRecoZ[2][i]) ); + + // Print the flash information for debugging + sFlashMatching += "Matching flash " + ProducerUtils::str(j) + " with time " + ProducerUtils::str(OpFlashTime[j]) + " and PE " + ProducerUtils::str(OpFlashPE[j]) + " in plane " + ProducerUtils::str(OpFlashPlane[j]) + " at distance " + ProducerUtils::str(OpFlashR) + " with residual " + ProducerUtils::str(OpFlashResidual) + "\n"; + + // Make a cut on the flash MaxPE/PE ratio and the number of hits + if ( OpFlashNHits[j] < fAdjOpFlashMinNHitCut || OpFlashMaxPE[j] / OpFlashPE[j] > fAdjOpFlashMaxPERatioCut ) { + continue; + } + + if ( lowe->SelectPDSFlashPE(TPCIDdriftTime[MVecTPC[2][i]], MVecTime[2][i] - OpFlashTime[j], MVecCharge[2][i], OpFlashPE[j]) ) { + // If the residual is smaller than the minimum residual, update the minimum residual and the matched flash. + if ( lowe->SelectPDSFlash(IsFirstFlash, TPCIDdriftTime[MVecTPC[2][i]], MVecTime[2][i], MVecCharge[2][i], MFlashTime, MFlashPE, OpFlashTime[j], OpFlashPE[j]) ) { + IsFirstFlash = false; + float a, b, c; + lowe->GetLightMapParameters("med", MVecCharge[2][i], a, b, c); + double RefPE = pow(10, a - a * b * (MVecTime[2][i] - OpFlashTime[j]) / TPCIDdriftTime[MVecTPC[2][i]] + c * pow((MVecTime[2][i] - OpFlashTime[j]) / TPCIDdriftTime[MVecTPC[2][i]], 2)); + MOpHitAmplitude = OpHitAmplitude[j]; + MFlashR = OpFlashR; + MFlashPE = OpFlashPE[j]; + MFlashFast = OpFlashFast[j]; + MFlashNHits = OpFlashNHits[j]; + MFlashPlane = OpFlashPlane[j]; + MFlashMaxPE = OpFlashMaxPE[j]; + MFlashPur = OpFlashPur[j]; + MFlashSTD = OpFlashSTD[j]; + MFlashTime = OpFlashTime[j]; + 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: " + 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]) + "\n" + + " - MainOpHitPE " + ProducerUtils::str(OpFlashMaxPE[j]) + " (PE); " + + " TotalPE " + ProducerUtils::str(OpFlashPE[j]) + " vs expected " + ProducerUtils::str(RefPE) + " (PE)\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" + + " - Found valid waveform: " + ProducerUtils::str(MFlashWaveformValid) + "\n"; + MatchedOpFlashX = MAdjFlashX; + // MatchedOpFlashResidual = OpFlashResidual; + MatchedOpFlashPE = MFlashPE; } } - OpFlashR = sqrt(pow(MAdjFlashX - OpFlashX[j], 2) + pow(MVecRecoY[2][i] - OpFlashY[j], 2)); + MAdjFlashResidual.push_back(OpFlashResidual); } - else if (fGeometry == "VD" && OpFlashPlane[j] == -1) { - sFlashMatching += "Skipping flash " + ProducerUtils::str(j) + " with unknown plane " + ProducerUtils::str(OpFlashPlane[j]) + "\n"; - continue; + + producer->PrintInColor(sFlashMatching, ProducerUtils::GetColor(sResultColor), "Debug"); + + if (fLabels[0] != "marley") { // If not marley, save the true signal particle information based on the mainID. + SignalParticlePDG = SignalPDGMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleE = SignalEMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleP = SignalPMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleK = SignalKMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleX = SignalStartXMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleY = SignalStartYMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleZ = SignalStartZMap[SignalPrimaryMap[MVecMainID[2][i]]]; + SignalParticleTime = SignalTimeMap[SignalPrimaryMap[MVecMainID[2][i]]]; } + sVertexReco += "*** Reconstructed Interaction Vertex: \n"; + sVertexReco += " - True X,Y,Z ( " + ProducerUtils::str(SignalStartXMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalStartYMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalStartZMap[SignalPrimaryMap[MVecMainID[2][i]]]) + " )" + "\n"; + sVertexReco += " - Main X,Y,Z ( " + ProducerUtils::str(SignalFinalXMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalFinalYMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalFinalZMap[SignalPrimaryMap[MVecMainID[2][i]]]) + " )" + "\n"; + sVertexReco += " - EDep X,Y,Z ( " + ProducerUtils::str(SignalMaxEDepXMap[MVecMainID[2][i]]) + ", " + ProducerUtils::str(SignalMaxEDepYMap[MVecMainID[2][i]]) + ", " + ProducerUtils::str(SignalMaxEDepZMap[MVecMainID[2][i]]) + " )" + "\n"; + sVertexReco += " - Reco X,Y,Z ( " + ProducerUtils::str(MatchedOpFlashX) + ", " + ProducerUtils::str(MVecRecoY[2][i]) + ", " + ProducerUtils::str(MVecRecoZ[2][i]) + " )"; + sClusterReco += sFlashReco; + sClusterReco += sVertexReco; - MAdjFlashR.push_back(OpFlashR); - MAdjFlashPE.push_back(OpFlashPE[j]); - MAdjFlashTime.push_back(OpFlashTime[j]); - MAdjFlashNHits.push_back(OpFlashNHits[j]); - MAdjFlashPlane.push_back(OpFlashPlane[j]); - MAdjFlashMaxPE.push_back(OpFlashMaxPE[j]); - MAdjFlashFast.push_back(OpFlashFast[j]); - MAdjFlashRecoX.push_back(OpFlashX[j]); - MAdjFlashRecoY.push_back(OpFlashY[j]); - MAdjFlashRecoZ.push_back(OpFlashZ[j]); - MAdjFlashSTD.push_back(OpFlashSTD[j]); - MAdjFlashPur.push_back(OpFlashPur[j]); - // Compute the residual between the predicted cluster signal and the flash - adjophits->FlashMatchResidual(OpFlashResidual, OpHitVec[j], MAdjFlashX, double(MVecRecoY[2][i]), double(MVecRecoZ[2][i])); - // Print the flash information for debugging - sFlashMatching += "Matching flash " + ProducerUtils::str(j) + " with time " + ProducerUtils::str(OpFlashTime[j]) + " and PE " + ProducerUtils::str(OpFlashPE[j]) + " in plane " + ProducerUtils::str(OpFlashPlane[j]) + " at distance " + ProducerUtils::str(OpFlashR) + " with residual " + ProducerUtils::str(OpFlashResidual) + "\n"; - // Make a cut on the flash MaxPE/PE ratio and the number of hits - if ( OpFlashNHits[j] < fAdjOpFlashMinNHitCut || OpFlashMaxPE[j] / OpFlashPE[j] > fAdjOpFlashMaxPERatioCut ) { + if (fClusterPreselectionFlashMatch && MatchedOpFlashPE < 0) { continue; - } - if ( lowe->SelectPDSFlashPE(TPCIDdriftTime[MVecTPC[2][i]], MVecCharge[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)) { - MFlashR = OpFlashR; - MFlashPE = OpFlashPE[j]; - MFlashFast = OpFlashFast[j]; - MFlashNHits = OpFlashNHits[j]; - MFlashPlane = OpFlashPlane[j]; - MFlashMaxPE = OpFlashMaxPE[j]; - MFlashPur = OpFlashPur[j]; - MFlashSTD = OpFlashSTD[j]; - MFlashTime = OpFlashTime[j]; - MFlashRecoX = OpFlashX[j]; - MFlashRecoY = OpFlashY[j]; - MFlashRecoZ = OpFlashZ[j]; - MFlashResidual = OpFlashResidual; - // Create an output string with the flash information. - sFlashReco = "*** Matched flash: \n - Purity " + ProducerUtils::str(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]) + - " Residual " + ProducerUtils::str(OpFlashResidual) + "\n" + - " - Reco Time,Y,Z ( " + ProducerUtils::str(MFlashTime) + ", " + ProducerUtils::str(OpFlashY[j]) + ", " + ProducerUtils::str(OpFlashZ[j]) + " )" + "\n"; - MatchedOpFlashX = MAdjFlashX; - MatchedOpFlashResidual = OpFlashResidual; - MatchedOpFlashPE = MFlashPE; - } - MAdjFlashResidual.push_back(OpFlashResidual); - } - - producer->PrintInColor(sFlashMatching, ProducerUtils::GetColor(sResultColor), "Debug"); - - if (fLabels[0] != "marley") { // If not marley, save the true signal particle information based on the mainID. - SignalParticlePDG = SignalPDGMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleE = SignalEMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleP = SignalPMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleK = SignalKMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleX = SignalStartXMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleY = SignalStartYMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleZ = SignalStartZMap[SignalPrimaryMap[MVecMainID[2][i]]]; - SignalParticleTime = SignalTimeMap[SignalPrimaryMap[MVecMainID[2][i]]]; - } - - sVertexReco += "*** Reconstructed Interaction Vertex: \n"; - sVertexReco += " - True X,Y,Z ( " + ProducerUtils::str(SignalStartXMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalStartYMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalStartZMap[SignalPrimaryMap[MVecMainID[2][i]]]) + " )" + "\n"; - sVertexReco += " - Main X,Y,Z ( " + ProducerUtils::str(SignalFinalXMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalFinalYMap[SignalPrimaryMap[MVecMainID[2][i]]]) + ", " + ProducerUtils::str(SignalFinalZMap[SignalPrimaryMap[MVecMainID[2][i]]]) + " )" + "\n"; - sVertexReco += " - EDep X,Y,Z ( " + ProducerUtils::str(SignalMaxEDepXMap[MVecMainID[2][i]]) + ", " + ProducerUtils::str(SignalMaxEDepYMap[MVecMainID[2][i]]) + ", " + ProducerUtils::str(SignalMaxEDepZMap[MVecMainID[2][i]]) + " )" + "\n"; - sVertexReco += " - Reco X,Y,Z ( " + ProducerUtils::str(MatchedOpFlashX) + ", " + ProducerUtils::str(MVecRecoY[2][i]) + ", " + ProducerUtils::str(MVecRecoZ[2][i]) + " )"; - sClusterReco += sFlashReco; - sClusterReco += sVertexReco; - - if (fClusterPreselectionFlashMatch && MatchedOpFlashPE < 0) { - continue; - } - // Fill the tree with the cluster information and the adjacent clusters and flashes. - MInd0Pur = MVecPur[0][i]; - MInd1Pur = MVecPur[1][i]; - MPur = MVecPur[2][i]; - MGen = MVecGen[i]; - MGenPur = MVecGenPur[i]; - MGenFrac = MVecGenFrac[i]; - MSignalFrac = {MVecFracE[i], MVecFracGa[i], MVecFracNe[i], MVecFracRest[i]}; - MInd0TPC = MVecTPC[0][i]; - MInd1TPC = MVecTPC[1][i]; - MTPC = MVecTPC[2][i]; - MInd0Charge = MVecCharge[0][i]; - MInd1Charge = MVecCharge[1][i]; - MCharge = MVecCharge[2][i]; - MInd0MaxCharge = MVecMaxCharge[0][i]; - MInd1MaxCharge = MVecMaxCharge[1][i]; - MMaxCharge = MVecMaxCharge[2][i]; - MInd0NHits = MVecNHits[0][i]; - MInd1NHits = MVecNHits[1][i]; - MNHit = MVecNHits[2][i]; - MInd0dTime = MVecdT[0][i]; - MInd1dTime = MVecdT[1][i]; - MTime = MVecTime[2][i]; - MRecX = MatchedOpFlashX; - MInd0RecoY = MVecRecoY[0][i]; - MInd1RecoY = MVecRecoY[1][i]; - MRecY = MVecRecoY[2][i]; - MRecZ = MVecRecoZ[2][i]; - MMainID = MVecMainID[2][i]; - // If mother exists add the mother information - const simb::MCParticle *MClTruth; - - int TerminalOutput = ProducerUtils::supress_stdout(); - MClTruth = pi_serv->TrackIdToParticle_P(MVecMainID[2][i]); - ProducerUtils::resume_stdout(TerminalOutput); - - if (MClTruth == 0) { - MMainVertex = {-1e6, -1e6, -1e6}; - MEndVertex = {-1e6, -1e6, -1e6}; - MMainPDG = 0; - MMainE = -1e6; - MMainP = -1e6; - MMainK = -1e6; - MMainTime = -1e6; - MMainParentVertex = {-1e6, -1e6, -1e6}; - MMainParentPDG = 0; - MMainParentE = -1e6; - MMainParentP = -1e6; - MMainParentK = -1e6; - MMainParentTime = -1e6; - } - else { - if (MFlashPur > 0) { - MFlashCorrect = true; - } - MMainVertex = {MClTruth->Vx(), MClTruth->Vy(), MClTruth->Vz()}; - MEndVertex = {MClTruth->EndX(), MClTruth->EndY(), MClTruth->EndZ()}; - MMainPDG = MClTruth->PdgCode(); - MMainE = 1e3 * MClTruth->E(); - MMainP = 1e3 * MClTruth->P(); - MMainK = MMainE - 1e3 * MClTruth->Mass(); - MMainTime = MClTruth->T(); - // If exists add the parent information. - const simb::MCParticle *MClParentTruth; - + } + // Fill the tree with the cluster information and the adjacent clusters and flashes. + MInd0Pur = MVecPur[0][i]; + MInd1Pur = MVecPur[1][i]; + MPur = MVecPur[2][i]; + MGen = MVecGen[i]; + MGenPur = MVecGenPur[i]; + MGenFrac = MVecGenFrac[i]; + MSignalFrac = {MVecFracE[i], MVecFracGa[i], MVecFracNe[i], MVecFracRest[i]}; + MInd0TPC = MVecTPC[0][i]; + MInd1TPC = MVecTPC[1][i]; + MTPC = MVecTPC[2][i]; + MInd0Charge = MVecCharge[0][i]; + MInd1Charge = MVecCharge[1][i]; + MCharge = MVecCharge[2][i]; + MInd0MaxCharge = MVecMaxCharge[0][i]; + MInd1MaxCharge = MVecMaxCharge[1][i]; + MMaxCharge = MVecMaxCharge[2][i]; + MInd0NHits = MVecNHits[0][i]; + MInd1NHits = MVecNHits[1][i]; + MNHit = MVecNHits[2][i]; + MInd0dTime = MVecdT[0][i]; + MInd1dTime = MVecdT[1][i]; + MTime = MVecTime[2][i]; + MRecX = MatchedOpFlashX; + MInd0RecoY = MVecRecoY[0][i]; + MInd1RecoY = MVecRecoY[1][i]; + MRecY = MVecRecoY[2][i]; + MRecZ = MVecRecoZ[2][i]; + MMainID = MVecMainID[2][i]; + const simb::MCParticle *MClTruth; // If mother exists add the mother information + int TerminalOutput = ProducerUtils::supress_stdout(); - MClParentTruth = pi_serv->TrackIdToParticle_P(MClTruth->Mother()); + MClTruth = pi_serv->TrackIdToParticle_P(MVecMainID[2][i]); ProducerUtils::resume_stdout(TerminalOutput); - if (MClParentTruth == 0) { + if (MClTruth == 0) { + MMainVertex = {-1e6, -1e6, -1e6}; + MEndVertex = {-1e6, -1e6, -1e6}; + MMainPDG = 0; + MMainE = -1e6; + MMainP = -1e6; + MMainK = -1e6; + MMainTime = -1e6; MMainParentVertex = {-1e6, -1e6, -1e6}; MMainParentPDG = 0; MMainParentE = -1e6; @@ -1836,41 +1921,70 @@ namespace solar MMainParentTime = -1e6; } else { - MMainParentVertex = {MClParentTruth->Vx(), MClParentTruth->Vy(), MClParentTruth->Vz()}; - MMainParentPDG = MClParentTruth->PdgCode(); - MMainParentE = 1e3 * MClParentTruth->E(); - MMainParentP = 1e3 * MClParentTruth->P(); - MMainParentK = MMainParentE - 1e3 * MClParentTruth->Mass(); - MMainParentTime = MClParentTruth->T(); + if ( MFlashPur > 0 ) { MFlashCorrect = true; } + MMainVertex = {MClTruth->Vx(), MClTruth->Vy(), MClTruth->Vz()}; + MEndVertex = {MClTruth->EndX(), MClTruth->EndY(), MClTruth->EndZ()}; + MMainPDG = MClTruth->PdgCode(); + MMainE = 1e3 * MClTruth->E(); + MMainP = 1e3 * MClTruth->P(); + MMainK = MMainE - 1e3 * MClTruth->Mass(); + MMainTime = MClTruth->T(); + const simb::MCParticle *MClParentTruth; // If exists add the parent information. + + int TerminalOutput = ProducerUtils::supress_stdout(); + MClParentTruth = pi_serv->TrackIdToParticle_P(MClTruth->Mother()); + ProducerUtils::resume_stdout(TerminalOutput); + + if (MClParentTruth == 0) { + MMainParentVertex = {-1e6, -1e6, -1e6}; + MMainParentPDG = 0; + MMainParentE = -1e6; + MMainParentP = -1e6; + MMainParentK = -1e6; + MMainParentTime = -1e6; + } + else { + MMainParentVertex = {MClParentTruth->Vx(), MClParentTruth->Vy(), MClParentTruth->Vz()}; + MMainParentPDG = MClParentTruth->PdgCode(); + MMainParentE = 1e3 * MClParentTruth->E(); + MMainParentP = 1e3 * MClParentTruth->P(); + MMainParentK = MMainParentE - 1e3 * MClParentTruth->Mass(); + MMainParentTime = MClParentTruth->T(); + } + } + + if (SignalParticleK < fMaxSignalK && SignalParticleK > 0) { + fSolarNuAnaTree->Fill(); + hDriftTime->Fill(MainElectronEndPointX, MTime); + hXTruth->Fill(MRecX - SignalParticleX, SignalParticleX); + hYTruth->Fill(MRecY - SignalParticleY, SignalParticleY); + hZTruth->Fill(MRecZ - SignalParticleZ, SignalParticleZ); } } - if (SignalParticleK < fMaxSignalK) { - fSolarNuAnaTree->Fill(); - hDriftTime->Fill(MainElectronEndPointX, MTime); - hXTruth->Fill(MRecX - SignalParticleX, SignalParticleX); - hYTruth->Fill(MRecY - SignalParticleY, SignalParticleY); - hZTruth->Fill(MRecZ - SignalParticleZ, SignalParticleZ); - } - } - // Check if the string sClusterReco is not empty and print it in color. - if (sClusterReco != "") { producer->PrintInColor(sClusterReco, ProducerUtils::GetColor(sResultColor)); } - } // Loop over clusters - if (SignalParticleK < fMaxSignalK) { + + // Check if the string sClusterReco is not empty and print it in color. + if (sClusterReco != "") { producer->PrintInColor(sClusterReco, ProducerUtils::GetColor(sResultColor)); } + } // Loop over clusters + fMCTruthTree->Fill(); SelectedEvents.push_back(1); - } + + + producer->PrintInColor("-----------------------------------------------------------------------------------------\n", ProducerUtils::GetColor("green")); + } // SelectedEvent + else { SelectedEvents.push_back(0); } - producer->PrintInColor("-----------------------------------------------------------------------------------------\n", ProducerUtils::GetColor("green")); - } // SelectedEvent } + void SolarNuAna::endJob() { + AnalyzedEvents = SelectedEvents.size(); producer->PrintInColor("Finished running the SolarNuAna module", ProducerUtils::GetColor("magenta")); fConfigTree->Fill(); } - //......................................................................................................................// + // Reset variables for each event void SolarNuAna::ResetVariables() { @@ -1892,19 +2006,24 @@ namespace solar MFlashRecoY = -1e6; MFlashRecoZ = -1e6; MFlashResidual = -1e6; + MFlashWaveform = {}; MFlashCorrect = false; - SignalParticleE = 0; - SignalParticleP = 0; - SignalParticleK = 0; - SignalParticleX = 0; - SignalParticleY = 0; - SignalParticleZ = 0; + SignalParticleE = -1e6; + SignalParticleP = -1e6; + SignalParticleK = -1e6; + SignalParticleX = -1e6; + SignalParticleY = -1e6; + SignalParticleZ = -1e6; SignalParticlePDG = 0; - SignalParticleTime = 0; + SignalParticleTime = -1e6; + 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..18d82707 100644 --- a/duneana/SolarNuAna/fcl/SolarNuAna.fcl +++ b/duneana/SolarNuAna/fcl/SolarNuAna.fcl @@ -1,3 +1,4 @@ +#include "LightMapData.fcl" #include "generator_labels.fcl" #include "SolarOpFlash.fcl" @@ -7,12 +8,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. @@ -24,14 +27,14 @@ BEGIN_PROLOG BackgroundLabelVector: [] #====================================================================================# - ClusterChargeVariable: "Integral" # Charge variable to use for cluster matching: "Integral" or "SummedADC". - OpHitTimeVariable: "PeakTime" # Time variable to use for OpHit selection: "PeakTime" or "StartTime". - XACathodeX: 0 # X position of the VD cathode XAs in [cm]. - XAMembraneY: 0 # Y position of the VD membrane XAs in [cm]. - XAFinalCapZ: 0 # Z position of the VD final cap XAs in [cm]. - XAStartCapZ: 0 # Z position of the VD start cap XAs in [cm]. + ClusterChargeVariable: "Integral" # Charge variable to use for cluster matching: "Integral" or "SummedADC". + OpHitTimeVariable: "PeakTime" # Time variable to use for OpHit selection: "PeakTime" or "StartTime". + XACathodeX: 0 # X position of the VD cathode XAs in [cm]. + XAMembraneY: 0 # Y position of the VD membrane XAs in [cm]. + XAFinalCapZ: 0 # Z position of the VD final cap XAs in [cm]. + XAStartCapZ: 0 # Z position of the VD start cap XAs in [cm]. ClusterAlgoTime: 12.5 # Time window to look for plane clusters in [us] units. - ClusterAlgoAdjChannel: 3 # Number of adjacent channels to look for plane clusters. + ClusterAlgoAdjChannel: 3 # Number of adjacent channels to look for plane clusters. GenerateSolarCluster: true # Generate SolarClusters from here defined matching cireteria ClusterMatchNHit: 2 # NHit fraction to match clusters. abs(NHitsCol - NHitsInd) / NHitsCol < ClusterMatchNHit. @@ -53,40 +56,44 @@ 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.16 # Negative time window to look for adj. OpHits in [us] units. + OpFlashAlgoMaxTime: 0.32 # 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. - 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. - AdjOpFlashMinPECut: 0.0 # Cut on the minimum OpHit PE. - AdjOpFlashMinPEAttenuation: 1.0 # Asymptotic growth attenuation on the minimum PE cut after a full drift length. - AdjOpFlashMinPEAttenuate: "flat" # Type of attenuation to apply to the minimum PE cut: "linear" or "asymptotic". - AdjOpFlashMinPEAttenuationStrength: 10 # For "asymptotic" attenuation. Strength expressed as powers of 10. - AdjOpFlashMinPELightMap: [ - ["Amplitude", [-1.1011e-7, 5.0670e-4, 1.6480e0]], - ["Attenuation", [1.7061e-7, -7.2303e-4, 1.0623e0]], - ["Correction", [3.1135e-7, -1.3050e-3, 7.6113e-1]] - ] # Quadratic fit parameters for attenuation correction based in logarithm light map. - FlashMatchByResidual: false # Match flashes by residual. Alternative is to match by MaxFlashPE. - + 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: 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: true # If true, the OpFlash matching is projected on the end cap planes for VD. + AdjOpFlashMaxPERatioCut: 1.0 # Cut on the maximum OpHit PE contribution to the total OpFlash PE. + AdjOpFlashMinPECut: 0.0 # Cut on the minimum OpHit PE. + AdjOpFlashMaxPECut: 1e9 # Cut on the maximum OpHit PE. + AdjOpFlashMinPEAttenuation: 0.0 # Asymptotic growth attenuation on the minimum PE cut after a full drift length. + AdjOpFlashMaxPEAttenuation: 0.0 # Asymptotic growth attenuation on the maximum PE cut after a full drift length. + AdjOpFlashMinPEAttenuate: "flat" # Type of attenuation to apply to the minimum PE cut: "linear" or "asymptotic". + AdjOpFlashMaxPEAttenuate: "flat" # Type of attenuation to apply to the maximum PE cut: "linear" or "asymptotic". + AdjOpFlashMinPEAttenuationStrength: 0 # For "asymptotic" attenuation. Strength expressed as powers of 10. + AdjOpFlashMaxPEAttenuationStrength: 0 # For "asymptotic" attenuation. Strength expressed as powers of 10. + AdjOpFlashMinPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashMinPELightMap + AdjOpFlashMaxPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashMaxPELightMap + AdjOpFlashPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashPELightMap + FlashMatchBy: "maximum" # Match flashes by residual. Alternative is to match by MaxFlashPE. + FlashMatchByPELightMapExponent: 2 # Exponent for PE light map matching. + SaveSignalDaughters: false # Save the daughters of the signal particles. SaveSignalEDep: false # Save the energy deposition of the marley energy deposition in the LAr. SaveSignalOpHits: false # Save the OpHits that are associated with the signal. SaveOpFlashInfo: false # Save the AdjOpFlash information for each preselection cluster. - SaveAdjOpFlashInfo: true # Save the AdjOpFlash information for each preselection cluster. + SaveAdjOpFlashInfo: true # Save the AdjOpFlash information for each preselection cluster. SaveTrackInfo: false # Save the MatchedTrack information for each preselection cluster. - TestLatestFeatures: true # If true, the latest features are used. + TestLatestFeatures: true # If true, the latest features are used. } legacy_solar_nu_ana_hd_v2: @local::solar_nu_ana_hd @@ -102,24 +109,23 @@ BEGIN_PROLOG solar_nu_ana_hd_lateralAPA.BackgroundLabelVector: @local::generator_dune10kt_1x2x6_lateralAPA solar_nu_ana_vd: @local::solar_nu_ana_hd_centralAPA - solar_nu_ana_vd.HitLabel: "gaushit" - solar_nu_ana_vd.OpHitLabel: "ophit10ppm" - solar_nu_ana_vd.OpFlashLabel: "solarflash" - solar_nu_ana_vd.OpFlashTimeOffset: 0 - solar_nu_ana_vd.BackgroundLabelVector: @local::generator_dunevd10kt_1x8x14_3view_30deg - solar_nu_ana_vd.OpHitTimeVariable: "StartTime" - solar_nu_ana_vd.XACathodeX: -327.5 - 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.OpFlashAlgoRad: 600. - solar_nu_ana_vd.AdjOpFlashMinPELightMap: [ - ["Amplitude", [1.7578e-7, -4.9589e-4, 1.2519e0]], - ["Attenuation", [4.9682e-7, -2.5242e-3, 1.2541e0]], - ["Correction", [1.0516e-7, -1.1637e-3, 1.5118e0]] - ] + 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.0 + solar_nu_ana_vd.BackgroundLabelVector: @local::generator_dunevd10kt_1x8x14_3view_30deg + solar_nu_ana_vd.OpHitTimeVariable: "StartTime" + solar_nu_ana_vd.XACathodeX: -327.5 + 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: 1.00 + solar_nu_ana_vd.OpFlashAlgoMaxTime: 1.60 + solar_nu_ana_vd.OpFlashAlgoRad: 600.0 + solar_nu_ana_vd.AdjOpFlashMinPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMinPELightMap + solar_nu_ana_vd.AdjOpFlashMaxPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMaxPELightMap + solar_nu_ana_vd.AdjOpFlashPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashPELightMap solar_nu_ana_vd_1x8x14_optimistic: @local::solar_nu_ana_vd solar_nu_ana_vd_1x8x14_optimistic.BackgroundLabelVector: @local::generator_dunevd10kt_1x8x14_3view_30deg_optimistic 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/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl index e94a7275..8ccffebb 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl @@ -3,4 +3,7 @@ #include "solar_ana_dune10kt_1x2x6.fcl" -physics.analyzers.solarnuana: @local::solar_nu_ana_hd_centralAPA \ No newline at end of file +physics.analyzers.solarnuana: @local::solar_nu_ana_hd_centralAPA +physics.analyzers.solarnuana.AdjOpFlashMinPEAttenuate: "light_map" +physics.analyzers.solarnuana.AdjOpFlashMaxPEAttenuate: "light_map" +physics.analyzers.solarnuana.FlashMatchBy: "light_map" \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl index 6a3a2636..b7a3353e 100644 --- a/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl @@ -3,4 +3,7 @@ #include "solar_ana_dune10kt_1x2x6.fcl" -physics.analyzers.solarnuana: @local::solar_nu_ana_hd_lateralAPA \ No newline at end of file +physics.analyzers.solarnuana: @local::solar_nu_ana_hd_lateralAPA +physics.analyzers.solarnuana.AdjOpFlashMinPEAttenuate: "light_map" +physics.analyzers.solarnuana.AdjOpFlashMaxPEAttenuate: "light_map" +physics.analyzers.solarnuana.FlashMatchBy: "maximum" diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl new file mode 100644 index 00000000..14a76e42 --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl @@ -0,0 +1,6 @@ +# Specific configuration for investigating signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD + +#include "solar_ana_single_flash_radiological_decay0_dune10kt_1x2x6_centralAPA.fcl" + +physics.analyzers.solarnuana.SaveOpFlashInfo: false \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl b/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl new file mode 100644 index 00000000..11609e7c --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefd/solar_ana_single_adjflash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl @@ -0,0 +1,6 @@ +# Specific configuration for investigating signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD + +#include "solar_ana_single_flash_radiological_decay0_dune10kt_1x2x6_lateralAPA.fcl" + +physics.analyzers.solarnuana.SaveOpFlashInfo: false \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl index 36272336..be512bf1 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -4,3 +4,4 @@ #include "solar_ana_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: true +physics.analyzers.solarnuana.AdjOpFlashEndCapProjection: false diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl index 4c5ab3d2..eae4390c 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_adjflash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -4,3 +4,4 @@ #include "solar_ana_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: true +physics.analyzers.solarnuana.AdjOpFlashEndCapProjection: false 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.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl index c169487d..cab53439 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -6,7 +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.AdjOpFlashEndCapProjection: false physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 6 physics.analyzers.solarnuana.AdjOpFlashX: 100. physics.analyzers.solarnuana.AdjOpFlashY: 100. 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..d7b66ea0 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,11 @@ 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.AdjOpFlashEndCapProjection: false +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 3885d71d..4e0bc208 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,8 +6,11 @@ physics.analyzers.solarnuana.GenerateAdjOpFlash: true physics.analyzers.solarnuana.SaveOpFlashInfo: true physics.analyzers.solarnuana.SaveAdjOpFlashInfo: true -physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: false -physics.analyzers.solarnuana.AdjOpFlashMinNHitCut: 6 +physics.analyzers.solarnuana.OpFlashAlgoPE: 1.5 +physics.analyzers.solarnuana.OpFlashAlgoTriggerPE: 1.5 +physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: true +physics.analyzers.solarnuana.AdjOpFlashEndCapProjection: false +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_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl index d98188e9..4198c970 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -4,3 +4,4 @@ #include "solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: true +physics.analyzers.solarnuana.AdjOpFlashEndCapProjection: false diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl index 1eea2106..ac0e9096 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_flash_yprojected_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -4,3 +4,4 @@ #include "solar_ana_flash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" physics.analyzers.solarnuana.AdjOpFlashMembraneProjection: true +physics.analyzers.solarnuana.AdjOpFlashEndCapProjection: false diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_exponent2_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_exponent2_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl new file mode 100644 index 00000000..b5abde7b --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_exponent2_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -0,0 +1,7 @@ +# Specific configuration for signal+bkg events with coustom flash generation +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD + +#include "solar_ana_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" + +physics.analyzers.solarnuana.SaveSignalDaughters: true +physics.analyzers.solarnuana.FlashMatchByPELightMapExponent: 2 diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_maximum_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_maximum_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl new file mode 100644 index 00000000..f22a7b35 --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_marley_adjflash_maximum_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -0,0 +1,7 @@ +# Specific configuration for signal+bkg events with coustom flash generation +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD + +#include "solar_ana_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" + +physics.analyzers.solarnuana.SaveSignalDaughters: true +physics.analyzers.solarnuana.FlashMatchBy: "maximum" \ No newline at end of file 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 diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl index fde3b94c..000f4877 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -4,4 +4,8 @@ #include "solar_ana_dunevd10kt_1x8x14_3view_30deg.fcl" +physics.analyzers.solarnuana.AdjOpFlashMinPEAttenuate: "light_map" +physics.analyzers.solarnuana.AdjOpFlashMaxPEAttenuate: "light_map" +physics.analyzers.solarnuana.FlashMatchBy: "maximum" + services.TFileService.fileName: "solar_ana_dunevd10kt_1x8x14_3view_30deg_nominal_hist.root" diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl index 5615e316..931a4e47 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_optimistic.fcl @@ -5,5 +5,8 @@ #include "solar_ana_dunevd10kt_1x8x14_3view_30deg.fcl" physics.analyzers.solarnuana: @local::solar_nu_ana_vd_1x8x14_optimistic +physics.analyzers.solarnuana.AdjOpFlashMinPEAttenuate: "light_map" +physics.analyzers.solarnuana.AdjOpFlashMaxPEAttenuate: "light_map" +physics.analyzers.solarnuana.FlashMatchBy: "light_map" services.TFileService.fileName: "solar_ana_dunevd10kt_1x8x14_3view_30deg_optimistic_hist.root" \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl index 035fad2f..88345579 100644 --- a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -5,5 +5,8 @@ #include "solar_ana_dunevd10kt_1x8x14_3view_30deg.fcl" physics.analyzers.solarnuana: @local::solar_nu_ana_vd_1x8x14_shielded +physics.analyzers.solarnuana.AdjOpFlashMinPEAttenuate: "light_map" +physics.analyzers.solarnuana.AdjOpFlashMaxPEAttenuate: "light_map" +physics.analyzers.solarnuana.FlashMatchBy: "light_map" services.TFileService.fileName: "solar_ana_dunevd10kt_1x8x14_3view_30deg_shielded_hist.root" \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl new file mode 100644 index 00000000..455b5a4d --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl @@ -0,0 +1,7 @@ +# Specific configuration for investigating signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD-VD + +#include "solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg.fcl" + +physics.analyzers.solarnuana.SignalLabel: "generator" +physics.analyzers.solarnuana.SaveOpFlashInfo: false \ No newline at end of file diff --git a/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl new file mode 100644 index 00000000..a9387e9f --- /dev/null +++ b/duneana/SolarNuAna/fcl/dunefdvd/solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl @@ -0,0 +1,7 @@ +# Specific configuration for investigating signal events. +# Run SolarNuAna on the output of a standard reco workflow for DUNE FD-VD + +#include "solar_ana_single_adjflash_radiological_decay0_dunevd10kt_1x8x14_3view_30deg_shielded.fcl" + +physics.analyzers.solarnuana.SignalLabel: "generator" +physics.analyzers.solarnuana.SaveOpFlashInfo: false \ No newline at end of file