Merge branch 'release/v10_14_00_01'

Author: nathanielerowe

CMakeLists.txt

@@ -15,7 +15,7 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.14.00)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.14.00.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 

sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc

@@ -195,6 +195,7 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p)
 
     produces< std::vector<int> >("MonPulses");
     produces< std::vector<int> >("MonPulseSizes"); 
+    produces<int>("pairsOverThreshold");
 }
 
 void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
@@ -278,6 +279,9 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
         auto sizesPtr = std::make_unique<std::vector<int>>();
         evt.put(std::move(flatPtr), "MonPulses");
         evt.put(std::move(sizesPtr), "MonPulseSizes");
+
+        auto pairFlag = std::make_unique<int>(-1);
+        evt.put(std::move(pairFlag), "pairsOverThreshold");
         return;
     }
 
@@ -623,30 +627,27 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
     std::vector<int> pulseSizes;
     MonPulsesFlat.clear();
     pulseSizes.clear();
-    int TotalFlash = pmtwvfmVec->size()/((int)fn_caenboards*PMTPerBoard); // pmtwvfmVec = waveHandle ???
+    int pmt_caenboards = (int)fn_caenboards-1;
+    int TotalFlash = pmtwvfmVec->size()/(pmt_caenboards*PMTPerBoard);
+
+    int numPairsOverThreshold = 0;
     for (int FlashCounter=0; FlashCounter<TotalFlash; FlashCounter++)
     {
       int WaveIndex = FlashCounter*PMTPerBoard;
       int WaveformSize = (*pmtwvfmVec)[WaveIndex].size();
+      int pairThisFlash = 0;
       std::vector<int> *MonPulse = new std::vector<int>(WaveformSize);
-      fTriggerService->ConstructMonPulse(*pmtwvfmVec, fmon_threshold, MonPulse, FlashCounter);
+      fTriggerService->ConstructMonPulse(*pmtwvfmVec, fmon_threshold, MonPulse, FlashCounter, &pairThisFlash);
       //MonPulsesAll.push_back(std::move(MonPulse));
       MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end());
       pulseSizes.push_back(MonPulse->size());
+      numPairsOverThreshold = numPairsOverThreshold + pairThisFlash;
       delete MonPulse;
     }
     // make ptrs
     auto flatPtr = std::make_unique<std::vector<int>>(std::move(MonPulsesFlat));
     auto sizesPtr = std::make_unique<std::vector<int>>(std::move(pulseSizes));
-
-    /*std::unique_ptr< std::vector< std::vector<int> > >  MonPulsesPtr(std::make_unique< std::vector< std::vector<int> > > ());
-    for (auto &pulse : MonPulsesAll) {
-      MonPulsesPtr->reserve(MonPulsesPtr->size() + pulse.size());
-      std::move(pulse.begin(), pulse.end(), std::back_inserter(*MonPulsesPtr));
-    }
-    // clean up the vector
-    for (unsigned i = 0; i < MonPulsesAll.size(); i++) MonPulsesAll[i] = std::vector<int>();
-    */ 
+    auto pairFlag = std::make_unique<int>(numPairsOverThreshold);
 
     board_frag_v.clear();
 
@@ -663,6 +664,7 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
 
     evt.put(std::move(flatPtr), "MonPulses");
     evt.put(std::move(sizesPtr), "MonPulseSizes");
+    evt.put(std::move(pairFlag), "pairsOverThreshold");
 }
 
 void sbndaq::SBNDPMTDecoder::get_fragments(artdaq::Fragment & frag, std::vector<std::vector<artdaq::Fragment>> & board_frag_v){

sbndcode/Decoders/PMT/run_pmtdecoder.fcl

@@ -11,7 +11,7 @@ services: {
     MonWidth: 12
     TotalCAENBoards: 8
     PMTPerBoard: 15
-    Baseline: 14250
+    Baseline: 14250 #Run 1: 14250, Run 2: 14257
     MC: false
   }
 }

sbndcode/JobConfigurations/base/detsim_drops.fcl

@@ -25,7 +25,8 @@ detsim_drops: [ @sequence::g4_drops
                 , "drop sim::SimEnergyDeposits_ionandscintout__*"
                 , "drop sim::SimEnergyDeposits_ionandscint_*_*"
                 , "drop sim::SimEnergyDeposits_largeant_*_*"
-				, "keep sim::SimEnergyDeposits_ionandscint_priorSCE_*"]
+				, "keep sim::SimEnergyDeposits_ionandscint_priorSCE_*"
+                , "drop raw::OpDetWaveforms_opdaq_*_*"]
 
 END_PROLOG
 

sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder.fcl

@@ -41,28 +41,27 @@ physics:
 {
   producers:
   {
-    opdecopmt:     @local::SBNDOpDeconvolutionPMT
-    ophitpmt:  @local::SBNDDecoOpHitFinderPMT
-    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0
-    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1
+    opdecopmt:     @local::SBNDOpDeconvolutionPMT_realisticMC
+    ophitpmt:  @local::SBNDDecoOpHitFinderPMT_realisticMC
+    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0_realisticMC
+    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1_realisticMC
 
     opdecoxarapuca:     @local::SBNDOpDeconvolutionXARAPUCA
     ophitxarapuca:  @local::SBNDDecoOpHitFinderXArapuca
     opflashdecotpc0arapuca:   @local::SBNDDecoSimpleFlashTPC0Arapuca
     opflashdecotpc1arapuca:   @local::SBNDDecoSimpleFlashTPC1Arapuca
-
   }
 
   reco: [
     opdecopmt,
     ophitpmt,
     opflashtpc0,
-    opflashtpc1,
+    opflashtpc1
 
-    opdecoxarapuca,
-    ophitxarapuca,
-    opflashdecotpc0arapuca,
-    opflashdecotpc1arapuca
+    #opdecoxarapuca,
+    #ophitxarapuca,
+    #opflashdecotpc0arapuca,
+    #opflashdecotpc1arapuca
   ]
 
 

sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl

@@ -62,6 +62,7 @@ SBNDDecoSimpleFlashTPC0_realisticMC.OpHitInputTime: "RiseTime"
 SBNDDecoSimpleFlashTPC0_realisticMC.UseT0Tool: true
 SBNDDecoSimpleFlashTPC0_realisticMC.ReadoutDelay: 0 //cable time delay in us
 SBNDDecoSimpleFlashTPC0_realisticMC.CorrectLightPropagation: true
+SBNDDecoSimpleFlashTPC0_realisticMC.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_MC2.root"
 
 #TPC1
 SBNDDecoSimpleFlashTPC1_realisticMC: @local::SBNDSimpleFlashTPC1
@@ -73,7 +74,7 @@ SBNDDecoSimpleFlashTPC1_realisticMC.OpHitInputTime: "RiseTime"
 SBNDDecoSimpleFlashTPC1_realisticMC.UseT0Tool: true
 SBNDDecoSimpleFlashTPC1_realisticMC.ReadoutDelay: 0 //cable time delay in us
 SBNDDecoSimpleFlashTPC1_realisticMC.CorrectLightPropagation: true
-SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data_v2.root"
+SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_MC2.root"
 
 
 ###### XA IDEAL MC ######

sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl

@@ -103,7 +103,5 @@ SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Method:2
 
 #####OpHit finder for PMT deconvolved waveforms#####
 SBNDDecoOpHitFinderPMT_realisticMC: @local::SBNDDecoOpHitFinderPMT_data
-####SPE area must be 1./DecoWaveformPrecision
-SBNDDecoOpHitFinderPMT_realisticMC.InputModule: "pmtpulseoscillation"
 
-END_PROLOG
+END_PROLOG

sbndcode/OpDetSim/TriggerEmulationService.cc

@@ -17,7 +17,7 @@ namespace calib {
     : fMonWidth(pset.get<int>("MonWidth", 12)),
       fTotalCAENBoards(pset.get<int>("TotalCAENBoards", 8)),
       PMTPerBoard(pset.get<int>("PMTPerBoard", 15)),
-      Baseline(pset.get<int>("Baseline", 14250)),
+      Baseline(pset.get<int>("Baseline")),
       fMC(pset.get<bool>("MC", true))
   {}
 
@@ -29,9 +29,11 @@ namespace calib {
       int MonThreshold,
       std::vector<int> *MonPulse,
       int FlashCounter,
-      int *numPairsOverThreshold
+      int *numPairsOverThreshold,
+      std::vector<int> PMT_Channels
   )
   {
+
       // Loop over the entries in our waveform vector
       // We care about getting the pairing correct
 
@@ -44,22 +46,37 @@ namespace calib {
           }
 
           std::map<int, const raw::OpDetWaveform*> channel_to_waveform;
-          for (const auto& wvf : fWaveforms)
+          for (const auto& wvf : fWaveforms) {
               channel_to_waveform[wvf.ChannelNumber()] = &wvf;
+          }
 
           std::vector<int> Pair2 = { 6,   8,  10,  12,  14,  16,  36,  38,  40,  60,  62,  66,  68, 70,  84,  86,  88,  90,  92,  94, 114, 116, 118, 138, 140, 144, 146, 148, 162, 164, 168, 170, 172, 192, 194, 196, 216, 218, 220, 222, 224, 226, 240, 242, 246, 248, 250, 270, 272, 274, 294, 296, 298, 300, 302, 304};
           std::vector<int> Pair1 = { 7,   9,  11,  13,  15,  17,  37,  39,  41,  61,  63,  67,  69, 71,  85,  87,  89,  91,  93,  95, 115, 117, 119, 139, 141, 145, 147, 149, 163, 165, 169, 171, 173, 193, 195, 197, 217, 219, 221, 223, 225, 227, 241, 243, 247, 249, 251, 271, 273, 275, 295, 297, 299, 301, 303, 305};
           std::vector<int> Unpaired = {65,  64, 143, 142, 167, 166, 245, 244};
           std::set<int> used_channels;
 
           // resize
-          int ReadoutSize = fWaveforms[0].size();
+          int ReadoutSize;
+          if (PMT_Channels.empty()) { std::cout<<"Warning: Please provide PMT channels list."<<std::endl; ReadoutSize = fWaveforms[Pair1[0]].size(); }
+          else ReadoutSize = fWaveforms[PMT_Channels[0]].size();
           MonPulse->assign(ReadoutSize, 0);
 
           for (size_t i = 0; i < Pair1.size(); ++i) {
               int ch1 = Pair1[i];
               int ch2 = Pair2[i];
 
+              // check that ch1 and ch2 show up in PMT_Channels
+              if (PMT_Channels.empty()) std::cout<<"Warning: Please provide PMT channels list to check if channels are PMT channels."<<std::endl;
+              else {
+                  if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch1) != PMT_Channels.end())) {
+                      throw cet::exception("TriggerEmulationService") << "Paired channel " << ch1 << " is not PMT channel. Check Pairs list. Skipping..." << std::endl;
+                      continue; 
+                  } if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch2) != PMT_Channels.end())) { 
+                      throw cet::exception("TriggerEmulationService") << "Paired channel " << ch2 << " is not PMT channel. Check Pairs list. Skipping..." << std::endl;
+                      continue; 
+                  }
+              }
+ 
               // skip if either waveform is missing
               if (channel_to_waveform.count(ch1) == 0 || channel_to_waveform.count(ch2) == 0) continue;
               // skip if already processed
@@ -82,6 +99,16 @@ namespace calib {
           }
 
           for (int ch : Unpaired) { // Unpaired channels
+
+              // check that ch1 and ch2 show up in PMT_Channels
+              if (PMT_Channels.empty()) std::cout<<"Warning: Please provide PMT channels list to check if channels are PMT channels."<<std::endl;
+              else {
+                  if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch) != PMT_Channels.end())) {
+                      throw cet::exception("TriggerEmulationService") << "Unpaired channel " << ch << " is not PMT channel. Check list. Skipping..." << std::endl;
+                      continue; 
+                  }
+              }
+ 
               if (used_channels.count(ch)) continue;
               if (channel_to_waveform.count(ch) == 0) continue;
 
@@ -143,6 +170,9 @@ namespace calib {
                   CAENChannel=CAENChannel+ChannelStep;
               } //loop over channels
           } //loop over boards
+
+          if (numPairsOverThreshold) *numPairsOverThreshold = *std::max_element(MonPulse->begin(), MonPulse->end());
+
       } // data  
   } // ConstructMonPulse
 

sbndcode/OpDetSim/TriggerEmulationService.h

@@ -37,7 +37,8 @@
              int MonThreshold, 
              std::vector<int> *MonPulse, 
              int FlashCounter,
-             int *numPairsOverThreshold = nullptr
+             int *numPairsOverThreshold = nullptr,
+             std::vector<int> PMT_Channels={}
            );
 
            int getTotalCAENBoards() const { return fTotalCAENBoards; } 

sbndcode/OpDetSim/digi_pmt_sbnd.fcl

@@ -23,8 +23,8 @@ sbnd_digipmt_alg:
   CableTime:             135        #time delay of the 30 m long readout cable in ns
 
   # Digitizer simulation
-  PMTADCDynamicRange:    14745      #in ADC values
-  PMTBaseline:           14745      #in ADC
+  PMTADCDynamicRange:    14250      #in ADC values
+  PMTBaseline:           14250      #in ADC
   PMTBaselineRMS:        2.6        #in ADC
 
   #Sample noise from data
@@ -35,13 +35,13 @@ sbnd_digipmt_alg:
   PMTDarkNoiseRate:      1000.0     #in Hz
 
   # Detection efficiencies
-  PMTCoatedVUVEff_tpc0:       0.0315      #PMT coated detection efficiency for direct (VUV) light
-  PMTCoatedVISEff_tpc0:       0.03493    #PMT coated detection efficiency for reflected (VIS) light
-  PMTUncoatedEff_tpc0:        0.03382   #PMT uncoated detection efficiency
+  PMTCoatedVUVEff_tpc0:       0.0392      #PMT coated detection efficiency for direct (VUV) light
+  PMTCoatedVISEff_tpc0:       0.026    #PMT coated detection efficiency for reflected (VIS) light
+  PMTUncoatedEff_tpc0:        0.0357   #PMT uncoated detection efficiency
 
-  PMTCoatedVUVEff_tpc1:       0.027      #PMT coated detection efficiency for direct (VUV) light
-  PMTCoatedVISEff_tpc1:       0.03493      #PMT coated detection efficiency for reflected (VIS) light
-  PMTUncoatedEff_tpc1:        0.03382      #PMT uncoated detection efficiency
+  PMTCoatedVUVEff_tpc1:       0.0392      #PMT coated detection efficiency for direct (VUV) light
+  PMTCoatedVISEff_tpc1:       0.026      #PMT coated detection efficiency for reflected (VIS) light
+  PMTUncoatedEff_tpc1:        0.0357      #PMT uncoated detection efficiency
 
   # Simulate gain fluctuations
   # (comment-out to skip gain fluctuations simulation)