Skip to content

Commit a25a48b

Browse files
npallatnpallat
committed
Fix: bug fix for correct MonPulse length, addition of number of pairs over threshold vector to events, updated Baseline to Run 2 Baseline, and checks that we are not considering any XArapuca channels
1 parent c895ade commit a25a48b

8 files changed

Lines changed: 123 additions & 77 deletions

‎sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc‎

Lines changed: 13 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -195,6 +195,7 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p)
195195

196196
produces< std::vector<int> >("MonPulses");
197197
produces< std::vector<int> >("MonPulseSizes");
198+
produces<int>("pairsOverThreshold");
198199
}
199200

200201
void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
@@ -278,6 +279,9 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
278279
auto sizesPtr = std::make_unique<std::vector<int>>();
279280
evt.put(std::move(flatPtr), "MonPulses");
280281
evt.put(std::move(sizesPtr), "MonPulseSizes");
282+
283+
auto pairFlag = std::make_unique<int>(-1);
284+
evt.put(std::move(pairFlag), "pairsOverThreshold");
281285
return;
282286
}
283287

@@ -623,30 +627,27 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
623627
std::vector<int> pulseSizes;
624628
MonPulsesFlat.clear();
625629
pulseSizes.clear();
626-
int TotalFlash = pmtwvfmVec->size()/((int)fn_caenboards*PMTPerBoard); // pmtwvfmVec = waveHandle ???
630+
int pmt_caenboards = (int)fn_caenboards-1;
631+
int TotalFlash = pmtwvfmVec->size()/(pmt_caenboards*PMTPerBoard);
632+
633+
int numPairsOverThreshold = 0;
627634
for (int FlashCounter=0; FlashCounter<TotalFlash; FlashCounter++)
628635
{
629636
int WaveIndex = FlashCounter*PMTPerBoard;
630637
int WaveformSize = (*pmtwvfmVec)[WaveIndex].size();
638+
int pairThisFlash = 0;
631639
std::vector<int> *MonPulse = new std::vector<int>(WaveformSize);
632-
fTriggerService->ConstructMonPulse(*pmtwvfmVec, fmon_threshold, MonPulse, FlashCounter);
640+
fTriggerService->ConstructMonPulse(*pmtwvfmVec, fmon_threshold, MonPulse, FlashCounter, &pairThisFlash);
633641
//MonPulsesAll.push_back(std::move(MonPulse));
634642
MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end());
635643
pulseSizes.push_back(MonPulse->size());
644+
numPairsOverThreshold = numPairsOverThreshold + pairThisFlash;
636645
delete MonPulse;
637646
}
638647
// make ptrs
639648
auto flatPtr = std::make_unique<std::vector<int>>(std::move(MonPulsesFlat));
640649
auto sizesPtr = std::make_unique<std::vector<int>>(std::move(pulseSizes));
641-
642-
/*std::unique_ptr< std::vector< std::vector<int> > > MonPulsesPtr(std::make_unique< std::vector< std::vector<int> > > ());
643-
for (auto &pulse : MonPulsesAll) {
644-
MonPulsesPtr->reserve(MonPulsesPtr->size() + pulse.size());
645-
std::move(pulse.begin(), pulse.end(), std::back_inserter(*MonPulsesPtr));
646-
}
647-
// clean up the vector
648-
for (unsigned i = 0; i < MonPulsesAll.size(); i++) MonPulsesAll[i] = std::vector<int>();
649-
*/
650+
auto pairFlag = std::make_unique<int>(numPairsOverThreshold);
650651

651652
board_frag_v.clear();
652653

@@ -663,6 +664,7 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt)
663664

664665
evt.put(std::move(flatPtr), "MonPulses");
665666
evt.put(std::move(sizesPtr), "MonPulseSizes");
667+
evt.put(std::move(pairFlag), "pairsOverThreshold");
666668
}
667669

668670
void sbndaq::SBNDPMTDecoder::get_fragments(artdaq::Fragment & frag, std::vector<std::vector<artdaq::Fragment>> & board_frag_v){

‎sbndcode/Decoders/PMT/run_pmtdecoder.fcl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ services: {
1111
MonWidth: 12
1212
TotalCAENBoards: 8
1313
PMTPerBoard: 15
14-
Baseline: 14250
14+
Baseline: 14257 #Run 1: 14250, Run 2: 14257
1515
MC: false
1616
}
1717
}

‎sbndcode/OpDetSim/TriggerEmulationService.cc‎

Lines changed: 34 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ namespace calib {
1717
: fMonWidth(pset.get<int>("MonWidth", 12)),
1818
fTotalCAENBoards(pset.get<int>("TotalCAENBoards", 8)),
1919
PMTPerBoard(pset.get<int>("PMTPerBoard", 15)),
20-
Baseline(pset.get<int>("Baseline", 14250)),
20+
Baseline(pset.get<int>("Baseline", 14257)),
2121
fMC(pset.get<bool>("MC", true))
2222
{}
2323

@@ -29,9 +29,11 @@ namespace calib {
2929
int MonThreshold,
3030
std::vector<int> *MonPulse,
3131
int FlashCounter,
32-
int *numPairsOverThreshold
32+
int *numPairsOverThreshold,
33+
std::vector<int> PMT_Channels
3334
)
3435
{
36+
3537
// Loop over the entries in our waveform vector
3638
// We care about getting the pairing correct
3739

@@ -44,22 +46,37 @@ namespace calib {
4446
}
4547

4648
std::map<int, const raw::OpDetWaveform*> channel_to_waveform;
47-
for (const auto& wvf : fWaveforms)
49+
for (const auto& wvf : fWaveforms) {
4850
channel_to_waveform[wvf.ChannelNumber()] = &wvf;
51+
}
4952

5053
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};
5154
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};
5255
std::vector<int> Unpaired = {65, 64, 143, 142, 167, 166, 245, 244};
5356
std::set<int> used_channels;
5457

5558
// resize
56-
int ReadoutSize = fWaveforms[0].size();
59+
int ReadoutSize;
60+
if (PMT_Channels.empty()) { std::cout<<"Warning: Please provide PMT channels list."<<std::endl; ReadoutSize = fWaveforms[Pair1[0]].size(); }
61+
else ReadoutSize = fWaveforms[PMT_Channels[0]].size();
5762
MonPulse->assign(ReadoutSize, 0);
5863

5964
for (size_t i = 0; i < Pair1.size(); ++i) {
6065
int ch1 = Pair1[i];
6166
int ch2 = Pair2[i];
6267

68+
// check that ch1 and ch2 show up in PMT_Channels
69+
if (PMT_Channels.empty()) std::cout<<"Warning: Please provide PMT channels list to check if channels are PMT channels."<<std::endl;
70+
else {
71+
if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch1) != PMT_Channels.end())) {
72+
std::cout<<"Paired channel "<<ch1<<" is not PMT channel. Check Pairs list. Skipping..."<<std::endl;
73+
continue;
74+
} if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch2) != PMT_Channels.end())) {
75+
std::cout<<"Paired channel "<<ch2<<" is not PMT channel. Check Pairs list. Skipping..."<<std::endl;
76+
continue;
77+
}
78+
}
79+
6380
// skip if either waveform is missing
6481
if (channel_to_waveform.count(ch1) == 0 || channel_to_waveform.count(ch2) == 0) continue;
6582
// skip if already processed
@@ -82,6 +99,16 @@ namespace calib {
8299
}
83100

84101
for (int ch : Unpaired) { // Unpaired channels
102+
103+
// check that ch1 and ch2 show up in PMT_Channels
104+
if (PMT_Channels.empty()) std::cout<<"Warning: Please provide PMT channels list to check if channels are PMT channels."<<std::endl;
105+
else {
106+
if (!(std::find(PMT_Channels.begin(), PMT_Channels.end(), ch) != PMT_Channels.end())) {
107+
std::cout<<"Unpaired channel "<<ch<<" is not PMT channel. Check list. Skipping..."<<std::endl;
108+
continue;
109+
}
110+
}
111+
85112
if (used_channels.count(ch)) continue;
86113
if (channel_to_waveform.count(ch) == 0) continue;
87114

@@ -143,6 +170,9 @@ namespace calib {
143170
CAENChannel=CAENChannel+ChannelStep;
144171
} //loop over channels
145172
} //loop over boards
173+
174+
if (numPairsOverThreshold) *numPairsOverThreshold = *std::max_element(MonPulse->begin(), MonPulse->end());
175+
146176
} // data
147177
} // ConstructMonPulse
148178

‎sbndcode/OpDetSim/TriggerEmulationService.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,8 @@
3737
int MonThreshold,
3838
std::vector<int> *MonPulse,
3939
int FlashCounter,
40-
int *numPairsOverThreshold = nullptr
40+
int *numPairsOverThreshold = nullptr,
41+
std::vector<int> PMT_Channels={}
4142
);
4243

4344
int getTotalCAENBoards() const { return fTotalCAENBoards; }

‎sbndcode/OpDetSim/opDetDigitizerSBND_module.cc‎

Lines changed: 70 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -169,13 +169,12 @@ namespace opdet {
169169
// Required functions.
170170
void produce(art::Event & e) override;
171171
std::vector<raw::OpDetWaveform> sliceWaveforms(std::vector<raw::OpDetWaveform> fWaveforms,
172-
int WaveIndex,
172+
std::vector<int> PMT_Channels,
173173
std::vector<int> *MonPulse,
174174
int PairMultiplicityThreshold,
175175
double tickPeriod,
176176
int ticksPerSlice,
177-
float PercentTicksBeforeCross,
178-
int PMTPerBoard);
177+
float PercentTicksBeforeCross);
179178
std::vector<std::vector<int>> sliceMonPulse(std::vector<int> *MonPulse,
180179
int PairMultiplicityThreshold,
181180
int ticksPerSlice,
@@ -314,6 +313,8 @@ namespace opdet {
314313
produces< std::vector< raw::OpDetWaveform > >();
315314
produces<bool>("triggerEmulation");
316315
produces<int>("pairsOverThreshold");
316+
produces< std::vector<int> >("pairsOverThresholdVec");
317+
produces<int>("numSlices");
317318
produces< std::vector< raw::OpDetWaveform > >("slicedWaveforms");
318319
produces< std::vector<int> >("MonPulses");
319320
produces< std::vector<int> >("MonPulseSizes");
@@ -373,65 +374,67 @@ namespace opdet {
373374
art::ServiceHandle<art::TFileService> tfs;
374375
art::ServiceHandle<calib::TriggerEmulationService> fTriggerService;
375376

376-
int PMTPerBoard = fTriggerService->getPMTPerBoard();
377-
int fTotalCAENBoards = fTriggerService->getTotalCAENBoards();
377+
std::vector<int> numPairsOverThreshold;
378378

379-
int TotalFlash = fWaveforms.size() / (fTotalCAENBoards * PMTPerBoard);
379+
// get PMT channels
380+
std::vector<int> PMT_Channels;
381+
for (const auto& wvf : fWaveforms) {
382+
auto ch = wvf.ChannelNumber();
383+
if (map.isPDType(ch, "pmt_uncoated") || map.isPDType(ch, "pmt_coated")) PMT_Channels.push_back(ch);
384+
}
380385

381-
int numPairsOverThreshold = 0;
386+
if (PMT_Channels.empty()) std::cout<<"Error: PMT has *NO* PMT channels in fWaveforms"<<std::endl;
387+
int WaveIndex = PMT_Channels[0];
388+
int WaveformSize = fWaveforms[WaveIndex].size();
382389

383-
// Loop through by flash -> compatible with ConstructMonPulse logic
384-
for (int FlashCounter = 0; FlashCounter < TotalFlash; ++FlashCounter) {
385-
int WaveIndex = FlashCounter*PMTPerBoard;
386-
int WaveformSize = fWaveforms[WaveIndex].size();
390+
std::vector<int>* MonPulse = new std::vector<int>(WaveformSize, 0);
387391

388-
std::vector<int>* MonPulse = new std::vector<int>(WaveformSize, 0);
389-
390-
int pairThisFlash = 0;
391-
// Send 3ms waveforms to ConstructMonPulse
392-
fTriggerService->ConstructMonPulse(fWaveforms, MonThreshold, MonPulse, FlashCounter, &pairThisFlash);
393-
numPairsOverThreshold = numPairsOverThreshold + pairThisFlash;
392+
int pairThisFlash = 0;
393+
// Send 3ms waveforms to ConstructMonPulse
394+
fTriggerService->ConstructMonPulse(fWaveforms, MonThreshold, MonPulse, 0, &pairThisFlash, PMT_Channels);
395+
numPairsOverThreshold.push_back(pairThisFlash);
394396

395-
double tickPeriod = sampling_rate(clockData);
397+
double tickPeriod = sampling_rate(clockData);
396398

397-
std::vector<raw::OpDetWaveform> SlicedWaveforms = sliceWaveforms(fWaveforms, WaveIndex, MonPulse, PairMultiplicityThreshold, tickPeriod, ticksPerSlice, PercentTicksBeforeCross, PMTPerBoard);
398-
std::vector<std::vector<int>> SlicedMonPulse = sliceMonPulse(MonPulse, PairMultiplicityThreshold, ticksPerSlice, PercentTicksBeforeCross);
399+
std::vector<raw::OpDetWaveform> SlicedWaveforms = sliceWaveforms(fWaveforms, PMT_Channels, MonPulse, PairMultiplicityThreshold, tickPeriod, ticksPerSlice, PercentTicksBeforeCross);
400+
std::vector<std::vector<int>> SlicedMonPulse = sliceMonPulse(MonPulse, PairMultiplicityThreshold, ticksPerSlice, PercentTicksBeforeCross);
399401

400-
SlicedWaveformsAll.push_back(std::move(SlicedWaveforms));
401-
MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end());
402-
pulseSizes.push_back(MonPulse->size());
402+
int numSlices = SlicedMonPulse.size();
403+
SlicedWaveformsAll.push_back(std::move(SlicedWaveforms));
404+
MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end());
405+
pulseSizes.push_back(MonPulse->size());
403406

404-
if (SaveNewPlots) {
405-
// Save histograms
406-
// Sliced waveforms
407-
for (size_t j; j < SlicedWaveformsAll.size(); ++j) {
408-
std::stringstream plotname2;
409-
plotname2 << "Sliced_waveforms_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << j;
410-
PlotWaveforms(SlicedWaveformsAll[j], plotname2.str());
411-
}
412-
// Long MonPulse
407+
if (SaveNewPlots) {
408+
// Save histograms
409+
// Sliced waveforms
410+
for (size_t j; j < SlicedWaveformsAll.size(); ++j) {
411+
std::stringstream plotname2;
412+
plotname2 << "Sliced_waveforms_" << e.id().event() << "_Mon_" << MonThreshold << "_slice" << j;
413+
PlotWaveforms(SlicedWaveformsAll[j], plotname2.str());
414+
}
415+
// Long MonPulse
416+
std::stringstream histname;
417+
histname << "Long_event_" << e.id().event() << "_Mon_" << MonThreshold;
418+
TH1D* MonHist = tfs->make<TH1D>(histname.str().c_str(), histname.str().c_str(),
419+
MonPulse->size(), 0.0, MonPulse->size() - 1);
420+
for (size_t i = 0; i < MonPulse->size(); i++) {
421+
MonHist->SetBinContent(i + 1, (*MonPulse)[i]);
422+
}
423+
// Sliced MonPulse
424+
for (size_t idx = 0; idx < SlicedMonPulse.size(); ++idx) {
425+
auto const& vec = SlicedMonPulse[idx];
413426
std::stringstream histname;
414-
histname << "Long_event_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter;
427+
histname << "Sliced_event_" << e.id().event() << "_Mon_" << MonThreshold << "_slice" << idx;
428+
415429
TH1D* MonHist = tfs->make<TH1D>(histname.str().c_str(), histname.str().c_str(),
416-
MonPulse->size(), 0.0, MonPulse->size() - 1);
417-
for (size_t i = 0; i < MonPulse->size(); i++) {
418-
MonHist->SetBinContent(i + 1, (*MonPulse)[i]);
419-
}
420-
// Sliced MonPulse
421-
for (size_t idx = 0; idx < SlicedMonPulse.size(); ++idx) {
422-
auto const& vec = SlicedMonPulse[idx];
423-
std::stringstream histname;
424-
histname << "Sliced_event_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << idx;
425-
426-
TH1D* MonHist = tfs->make<TH1D>(histname.str().c_str(), histname.str().c_str(),
427-
vec.size(), 0.0, vec.size() - 1);
428-
for (size_t i = 0; i < vec.size(); i++) {
429-
MonHist->SetBinContent(i + 1, vec[i]);
430-
}
430+
vec.size(), 0.0, vec.size() - 1);
431+
for (size_t i = 0; i < vec.size(); i++) {
432+
MonHist->SetBinContent(i + 1, vec[i]);
431433
}
432-
}
434+
}
435+
}
433436
delete MonPulse;
434-
}
437+
//}
435438

436439
// find the trigger locations for the waveforms - old version, keeping for validation
437440
for (const raw::OpDetWaveform &waveform : fWaveforms) {
@@ -485,9 +488,21 @@ namespace opdet {
485488
e.put(std::move(triggerFlag), "triggerEmulation");
486489

487490

491+
// put number of slices in the event
492+
auto numSlicesFlag = std::make_unique<int>(numSlices);
493+
e.put(std::move(numSlicesFlag), "numSlices");
494+
495+
488496
// put trigger pair result in the event
489-
auto pairFlag = std::make_unique<int>(numPairsOverThreshold);
497+
int max_numPairsOverThreshold = 0;
498+
if (!numPairsOverThreshold.empty()) {
499+
for (int n : numPairsOverThreshold) if (n > max_numPairsOverThreshold) max_numPairsOverThreshold = n;
500+
}
501+
auto pairFlag = std::make_unique<int>(max_numPairsOverThreshold);
490502
e.put(std::move(pairFlag), "pairsOverThreshold");
503+
// put trigger pair result in the event
504+
auto pairFlagVec = std::make_unique<std::vector<int>>(numPairsOverThreshold);
505+
e.put(std::move(pairFlagVec), "pairsOverThresholdVec");
491506

492507

493508
// put sliced waveforms in the event
@@ -504,7 +519,6 @@ namespace opdet {
504519
// put the waveforms in the event
505520
e.put(std::move(SlicedWaveformsPtr), "slicedWaveforms");
506521

507-
508522
// put MonPulses in the event
509523
auto flatPtr = std::make_unique<std::vector<int>>(std::move(MonPulsesFlat));
510524
e.put(std::move(flatPtr), "MonPulses");
@@ -611,13 +625,12 @@ namespace opdet {
611625
// sliceWaveforms function
612626
std::vector<raw::OpDetWaveform> opDetDigitizerSBND::sliceWaveforms(
613627
std::vector<raw::OpDetWaveform> fWaveforms,
614-
int WaveIndex,
628+
std::vector<int> PMT_Channels,
615629
std::vector<int>* MonPulse,
616630
int PairMultiplicityThreshold,
617631
double tickPeriod,
618632
int ticksPerSlice,
619-
float PercentTicksBeforeCross,
620-
int PMTPerBoard)
633+
float PercentTicksBeforeCross)
621634
{
622635
// before and after crossing point (default is ~20% and ~80%)
623636
int ticksBeforeCross = static_cast<int>(std::round(PercentTicksBeforeCross*ticksPerSlice));
@@ -628,8 +641,8 @@ namespace opdet {
628641

629642
std::vector<raw::OpDetWaveform> SlicedWaveforms;
630643
// loop through channels
631-
for (int chan = 0; chan < PMTPerBoard; ++chan) {
632-
const raw::OpDetWaveform& wf = fWaveforms[WaveIndex + chan];
644+
for (int chan : PMT_Channels) {
645+
const raw::OpDetWaveform& wf = fWaveforms[chan];
633646

634647
for (auto [start, end] : intervals) {
635648
double sliceTime = wf.TimeStamp() + start * tickPeriod;

‎sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ sbnd_opdetdigitizer:
1515
ticksPerSlice: 5000 # corresponds to 10us
1616
PercentTicksBeforeCross: 0.2
1717
MonThreshold: 15
18-
PairMultiplicityThreshold: 15
18+
PairMultiplicityThreshold: 4
1919
SaveNewPlots: true
2020
SaveOldPlots: false
2121

‎sbndcode/OpDetSim/run_BeamRatesCalib.fcl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ services:
1717
MonWidth: 12
1818
TotalCAENBoards: 8
1919
PMTPerBoard: 15
20-
Baseline: 14250
20+
Baseline: 14257 #Run 1: 14250, Run 2: 14257
2121
MC: false
2222
}
2323
}

‎sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ sbnd_triggeremulation_service: {
66
MonWidth: 12
77
TotalCAENBoards: 8
88
PMTPerBoard: 15
9-
Baseline: 14250
9+
Baseline: 14257 #Run 1: 14250, Run 2: 14257
1010
MC: true
1111
}
1212

0 commit comments

Comments
 (0)