Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:20:24

0001 // use #include "" only for your local include and put
0002 // those in the first line(s) before any #include <>
0003 // otherwise you are asking for weird behavior
0004 // (more info - check the difference in include path search when using "" versus <>)
0005 
0006 #include "HcalMon.h"
0007 
0008 #include <onlmon/OnlMon.h>  // for OnlMon
0009 #include <onlmon/OnlMonDB.h>
0010 #include <onlmon/OnlMonServer.h>
0011 #include <onlmon/pseudoRunningMean.h>
0012 #include <onlmon/triggerEnum.h>
0013 
0014 #include <calobase/TowerInfoDefs.h>
0015 #include <caloreco/CaloWaveformFitting.h>
0016 
0017 #include <Event/Event.h>
0018 #include <Event/eventReceiverClient.h>
0019 #include <Event/msg_profile.h>
0020 
0021 #include <TH1.h>
0022 #include <TH2.h>
0023 #include <TProfile.h>
0024 #include <TProfile2D.h>
0025 
0026 #include <cmath>
0027 #include <cstdio>  // for printf
0028 #include <fstream>
0029 #include <iostream>
0030 #include <sstream>
0031 #include <string>  // for allocator, string, char_traits
0032 
0033 enum
0034 {
0035   TRGMESSAGE = 1,
0036   FILLMESSAGE = 2
0037 };
0038 
0039 HcalMon::HcalMon(const std::string& name)
0040   : OnlMon(name)
0041 {
0042   // leave ctor fairly empty, its hard to debug if code crashes already
0043   // during a new HcalMon()
0044   // if name start with O then packetlow = 8001, packethigh = 8008
0045   // if name start with I then packetlow = 7001, packethigh = 7008
0046   if (name[0] == 'O')
0047   {
0048     packetlow = 8001;
0049     packethigh = 8008;
0050   }
0051   else if (name[0] == 'I')
0052   {
0053     packetlow = 7001;
0054     packethigh = 7008;
0055   }
0056   else
0057   {
0058     std::cout << "HcalMon::HcalMon - unknown name(need to be OHCALMON or IHCALMON to know what packet to run) " << name << std::endl;
0059     exit(1);
0060   }
0061   return;
0062 }
0063 
0064 HcalMon::~HcalMon()
0065 {
0066   // you can delete NULL pointers it results in a NOOP (No Operation)
0067   delete WaveformProcessing;
0068   for (auto iter : rm_vector_sectAvg)
0069   {
0070     delete iter;
0071   }
0072   for (auto iter : rm_vector_twr)
0073   {
0074     delete iter;
0075   }
0076   for (auto iter : rm_packet_number)
0077   {
0078     delete iter;
0079   }
0080   for (auto iter : rm_packet_length)
0081   {
0082     delete iter;
0083   }
0084   for (auto iter : rm_packet_chans)
0085   {
0086     delete iter;
0087   }
0088   for (auto iter : rm_vector_twrTime)
0089   {
0090     delete iter;
0091   }
0092   for (auto iter : rm_vector_twrhit)
0093   {
0094     delete iter;
0095   }
0096   for (auto iter : rm_vector_twrhit_alltrig)
0097   {
0098     delete iter;
0099   }
0100 
0101   if (erc)
0102   {
0103     delete erc;
0104   }
0105 
0106   return;
0107 }
0108 
0109 const int depth = 200000;
0110 const int packet_depth = 1000;
0111 const int historyLength = 100;
0112 const int historyScaleDown = 100;
0113 // const int n_channel = 48;
0114 float hit_threshold = 30;
0115 float waveform_hit_threshold = 100;
0116 const int n_samples_show = 31;
0117 
0118 int HcalMon::Init()
0119 {
0120   // read our calibrations from HcalMonData.dat
0121   /*
0122 
0123   const char *hcalcalib = getenv("HCALCALIB");
0124   if (!hcalcalib)
0125   {
0126     std::cout << "HCALCALIB environment variable not set" << std::endl;
0127     exit(1);
0128   }
0129   std::string fullfile = std::string(hcalcalib) + "/" + "HcalMonData.dat";
0130   std::ifstream calib(fullfile);
0131   calib.close();
0132   */
0133   // use printf for stuff which should go the screen but not into the message
0134   // system (all couts are redirected)
0135   printf("doing the Init\n");
0136 
0137   for (int itrig = 0; itrig < 64; itrig++)
0138   {
0139     h2_hcal_hits_trig[itrig] = new TH2F(Form("h2_hcal_hits_trig_%d", itrig), "", 24, 0, 24, 64, 0, 64);
0140   }
0141 
0142   pr_zsFrac_etaphi = new TProfile2D("pr_zsFrac_etaphi", "", 24, 0, 24, 64, 0, 64);
0143   pr_zsFrac_etaphi_all = new TProfile2D("pr_zsFrac_etaphi_all", "", 24, 0, 24, 64, 0, 64);
0144   h_hcal_trig = new TH1F("h_hcal_trig", "", 64, 0, 64);
0145   h2_hcal_rm = new TH2F("h2_hcal_rm", "", 24, 0, 24, 64, 0, 64);
0146   h2_hcal_rm_alltrig = new TH2F("h2_hcal_rm_alltrig", "", 24, 0, 24, 64, 0, 64);
0147   h2_hcal_mean = new TH2F("h2_hcal_mean", "", 24, 0, 24, 64, 0, 64);
0148   h2_hcal_time = new TH2F("h2_hcal_time", "", 24, 0, 24, 64, 0, 64);
0149   h2_hcal_hits = new TH2F("h2_hcal_hits", "", 24, 0, 24, 64, 0, 64);
0150   h2_hcal_waveform = new TH2F("h2_hcal_waveform", "", n_samples_show, 0.5, n_samples_show + 0.5, 1000, 0, 15000);
0151   h2_hcal_correlation = new TH2F("h2_hcal_correlation", "", 200, 0, 100000, 200, 0, 150000);
0152   h_event = new TH1F("h_event", "", 1, 0, 1);
0153   h_waveform_twrAvg = new TH1F("h_waveform_twrAvg", "", n_samples_show, 0.5, n_samples_show + 0.5);
0154   h_waveform_time = new TH1F("h_waveform_time", "", n_samples_show, 0.5, n_samples_show + 0.5);
0155   h_waveform_pedestal = new TH1F("h_waveform_pedestal", "", 5e3, 0, 5e3);
0156   h_sectorAvg_total = new TH1F("h_sectorAvg_total", "", 32, 0.5, 32.5);
0157   // number of towers above threshold per event
0158   h_ntower = new TH1F("h_ntower", "", 100, 0, 800);
0159   // packet stuff
0160   h1_packet_number = new TH1F("h1_packet_number", "", 8, packetlow - 0.5, packethigh + 0.5);
0161   h1_packet_length = new TH1F("h1_packet_length", "", 8, packetlow - 0.5, packethigh + 0.5);
0162   h1_packet_chans = new TH1F("h1_packet_chans", "", 8, packetlow - 0.5, packethigh + 0.5);
0163   h1_packet_event = new TH1F("h1_packet_event", "", 8, packetlow - 0.5, packethigh + 0.5);
0164   h_caloPack_gl1_clock_diff = new TH2F("h_caloPack_gl1_clock_diff", "", 8, packetlow - 0.5, packethigh + 0.5, 65536, 0, 65536);
0165   h_evtRec = new TProfile("h_evtRec", "", 1, 0, 1);
0166 
0167   p2_pre_post = new TProfile2D("p2_pre_post", "", 24, 0, 24, 64, 0, 64, "S");
0168 
0169   for (int ih = 0; ih < Nsector; ih++)
0170   {
0171     h_rm_sectorAvg[ih] = new TH1F(Form("h_rm_sectorAvg_s%d", ih), "", historyLength, 0, historyLength * historyScaleDown);
0172   }
0173   for (int ieta = 0; ieta < 24; ieta++)
0174   {
0175     for (int iphi = 0; iphi < 64; iphi++)
0176     {
0177       h_rm_tower[ieta][iphi] = new TH1F(Form("h_rm_tower_%d_%d", ieta, iphi), Form("multiplicity running mean of tower ieta=%d, iphi=%d", ieta, iphi), historyLength, 0, historyLength * historyScaleDown);
0178     }
0179   }
0180   // make the per-packet running mean objects
0181   // 32 packets and 48 channels for hcal detectors
0182   for (int i = 0; i < Nsector; i++)
0183   {
0184     rm_vector_sectAvg.push_back(new pseudoRunningMean(1, depth));
0185   }
0186   for (int i = 0; i < Ntower; i++)
0187   {
0188     rm_vector_twr.push_back(new pseudoRunningMean(1, depth));
0189     rm_vector_twrTime.push_back(new pseudoRunningMean(1, depth));
0190     rm_vector_twrhit.push_back(new pseudoRunningMean(1, depth));
0191     rm_vector_twrhit_alltrig.push_back(new pseudoRunningMean(1, depth));
0192   }
0193   for (int i = 0; i < 8; i++)
0194   {
0195     rm_packet_number.push_back(new pseudoRunningMean(1, packet_depth));
0196     rm_packet_length.push_back(new pseudoRunningMean(1, packet_depth));
0197     rm_packet_chans.push_back(new pseudoRunningMean(1, packet_depth));
0198   }
0199 
0200   OnlMonServer* se = OnlMonServer::instance();
0201   // register histograms with server otherwise client won't get them
0202   se->registerHisto(this, h2_hcal_hits);
0203   for (int itrig = 0; itrig < 64; itrig++)
0204   {
0205     se->registerHisto(this, h2_hcal_hits_trig[itrig]);
0206   }
0207 
0208   se->registerHisto(this, pr_zsFrac_etaphi);
0209   se->registerHisto(this, pr_zsFrac_etaphi_all);
0210   se->registerHisto(this, h_hcal_trig);
0211   se->registerHisto(this, h_evtRec);
0212   se->registerHisto(this, h_caloPack_gl1_clock_diff);
0213   se->registerHisto(this, h2_hcal_rm);
0214   se->registerHisto(this, h2_hcal_rm_alltrig);
0215   se->registerHisto(this, h2_hcal_mean);
0216   se->registerHisto(this, h2_hcal_time);
0217   se->registerHisto(this, h2_hcal_waveform);
0218   se->registerHisto(this, h_event);
0219   se->registerHisto(this, h_sectorAvg_total);
0220   se->registerHisto(this, h_waveform_twrAvg);
0221   se->registerHisto(this, h_waveform_time);
0222   se->registerHisto(this, h_waveform_pedestal);
0223   se->registerHisto(this, h_ntower);
0224   se->registerHisto(this, h1_packet_number);
0225   se->registerHisto(this, h1_packet_length);
0226   se->registerHisto(this, h1_packet_chans);
0227   se->registerHisto(this, h1_packet_event);
0228   se->registerHisto(this, h2_hcal_correlation);
0229   se->registerHisto(this, p2_pre_post);
0230 
0231   for (auto& ih : h_rm_sectorAvg)
0232   {
0233     se->registerHisto(this, ih);
0234   }
0235 
0236   for (auto& ieta : h_rm_tower)
0237   {
0238     for (auto& iphi : ieta)
0239     {
0240       se->registerHisto(this, iphi);
0241     }
0242   }
0243 
0244   Reset();
0245 
0246   // initialize waveform extraction tool
0247   WaveformProcessing = new CaloWaveformFitting();
0248 
0249   std::string hcaltemplate;
0250   if (getenv("HCALCALIB"))
0251   {
0252     hcaltemplate = getenv("HCALCALIB");
0253   }
0254   else
0255   {
0256     hcaltemplate = ".";
0257   }
0258   hcaltemplate += std::string("/testbeam_ohcal_template.root");
0259   // WaveformProcessing->initialize_processing(hcaltemplate);
0260 
0261   if (anaGL1)
0262   {
0263     erc = new eventReceiverClient("gl1daq");
0264   }
0265 
0266   return 0;
0267 }
0268 
0269 std::vector<float> HcalMon::getSignal(Packet* p, const int channel)
0270 {
0271   double baseline = 0;
0272   for (int s = 0; s < 3; s++)
0273   {
0274     baseline += p->iValue(s, channel);
0275   }
0276   baseline /= 3.;
0277 
0278   double signal = 0;
0279   int sample = 0;
0280   for (int s = 3; s < p->iValue(0, "SAMPLES"); s++)
0281   {
0282     if (signal > p->iValue(s, channel))
0283     {
0284       signal = p->iValue(s, channel);
0285       sample = s;
0286     }
0287   }
0288   signal -= baseline;
0289 
0290   std::vector<float> result;
0291   result.push_back(signal);
0292   result.push_back(sample);
0293   result.push_back(baseline);
0294 
0295   return result;
0296 }
0297 
0298 std::vector<float> HcalMon::anaWaveform(Packet* p, const int channel)
0299 {
0300   std::vector<float> waveform;
0301   // waveform.reserve(p->iValue(0, "SAMPLES"));
0302   float supppressed = 1;
0303   if (p->iValue(channel, "SUPPRESSED"))
0304   {
0305     waveform.push_back(p->iValue(channel, "PRE"));
0306     waveform.push_back(p->iValue(channel, "POST"));
0307   }
0308   else
0309   {
0310     supppressed = 0;
0311     for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0312     {
0313       waveform.push_back(p->iValue(s, channel));
0314     }
0315   }
0316   std::vector<std::vector<float>> multiple_wfs;
0317   multiple_wfs.push_back(waveform);
0318 
0319   std::vector<std::vector<float>> fitresults_ohcal;
0320   // fitresults_ohcal = WaveformProcessing->process_waveform(multiple_wfs);
0321   fitresults_ohcal = WaveformProcessing->calo_processing_fast(multiple_wfs);
0322 
0323   std::vector<float> result;
0324   result = fitresults_ohcal.at(0);
0325   result.push_back(supppressed);
0326 
0327   return result;
0328 }
0329 
0330 int HcalMon::BeginRun(const int /* runno */)
0331 {
0332   // reset the thresholds
0333   hit_threshold = 30;
0334   waveform_hit_threshold = 100;
0335 
0336   // if you need to read calibrations on a run by run basis
0337   // this is the place to do it
0338 
0339   std::vector<runningMean*>::iterator rm_it;
0340   for (rm_it = rm_vector_sectAvg.begin(); rm_it != rm_vector_sectAvg.end(); ++rm_it)
0341   {
0342     (*rm_it)->Reset();
0343   }
0344   for (rm_it = rm_vector_twr.begin(); rm_it != rm_vector_twr.end(); ++rm_it)
0345   {
0346     (*rm_it)->Reset();
0347   }
0348   for (rm_it = rm_packet_number.begin(); rm_it != rm_packet_number.end(); ++rm_it)
0349   {
0350     (*rm_it)->Reset();
0351   }
0352   for (rm_it = rm_packet_length.begin(); rm_it != rm_packet_length.end(); ++rm_it)
0353   {
0354     (*rm_it)->Reset();
0355   }
0356   for (rm_it = rm_packet_chans.begin(); rm_it != rm_packet_chans.end(); ++rm_it)
0357   {
0358     (*rm_it)->Reset();
0359   }
0360   for (rm_it = rm_vector_twrTime.begin(); rm_it != rm_vector_twrTime.end(); ++rm_it)
0361   {
0362     (*rm_it)->Reset();
0363   }
0364   for (rm_it = rm_vector_twrhit.begin(); rm_it != rm_vector_twrhit.end(); ++rm_it)
0365   {
0366     (*rm_it)->Reset();
0367   }
0368   for (rm_it = rm_vector_twrhit_alltrig.begin(); rm_it != rm_vector_twrhit_alltrig.end(); ++rm_it)
0369   {
0370     (*rm_it)->Reset();
0371   }
0372   if (anaGL1)
0373   {
0374     OnlMonServer* se = OnlMonServer::instance();
0375     se->UseGl1();
0376   }
0377   return 0;
0378 }
0379 
0380 int HcalMon::process_event(Event* e /* evt */)
0381 {
0382   if (e->getEvtType() >= 8)  /// special event where we do not read out the calorimeters
0383   {
0384     return 0;
0385   }
0386   evtcnt++;
0387   h_waveform_twrAvg->Reset();  // only record the latest event waveform
0388   h1_packet_event->Reset();
0389   unsigned int towerNumber = 0;
0390   float sectorAvg[Nsector] = {0};
0391   int npacket1 = 0;
0392   int npacket2 = 0;
0393   float energy1 = 0;
0394   float energy2 = 0;
0395 
0396   bool fillhist = true;
0397   std::vector<bool> trig_bools;
0398   trig_bools.resize(64);
0399   long long int gl1_clock = 0;
0400   bool have_gl1 = false;
0401   if (anaGL1)
0402   {
0403     int evtnr = e->getEvtSequence();
0404     Event* gl1Event = erc->getEvent(evtnr);
0405     if (gl1Event)
0406     {
0407       OnlMonServer* se = OnlMonServer::instance();
0408       se->IncrementGl1FoundCounter();
0409       have_gl1 = true;
0410       Packet* p = gl1Event->getPacket(14001);
0411       h_evtRec->Fill(0.0, 1.0);
0412       if (p)
0413       {
0414         gl1_clock = p->lValue(0, "BCO");
0415         uint64_t triggervec = p->lValue(0, "ScaledVector");
0416         for (int i = 0; i < 64; i++)
0417         {
0418           bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0419           trig_bools[i] = trig_decision;
0420 
0421           if (trig_decision)
0422           {
0423             h_hcal_trig->Fill(i);
0424           }
0425           triggervec = (triggervec >> 1U);
0426         }
0427 
0428         delete p;
0429       }
0430       delete gl1Event;
0431     }
0432     else
0433     {
0434       std::cout << "GL1 event is null" << std::endl;
0435       h_evtRec->Fill(0.0, 0.0);
0436     }
0437     // this is for only fill certain histogram with the MBD>=1 trigger or cosmic single trigger(and they should never run together!!)
0438     if (usetrig4_10)
0439     {
0440       // commenting out until we have a better way to handle cosmic bits -- tanner
0441       if (trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX10) == 0 && trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX30) == 0 && trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX150) == 0 && trig_bools.at(TriggerEnum::BitCodes::HCAL_SINGLES) == 0)
0442       {
0443         fillhist = false;
0444       }
0445       // if we have hcal single cosmic trigger we are in cosmic running mode and need to adjust thresholds accordingly
0446       if (trig_bools.at(TriggerEnum::BitCodes::RANDOM) || trig_bools.at(TriggerEnum::BitCodes::HCAL_SINGLES) || trig_bools.at(TriggerEnum::BitCodes::HCAL_NARROW_VERT) || trig_bools.at(TriggerEnum::BitCodes::HCAL_WIDE_VERT) || trig_bools.at(TriggerEnum::BitCodes::HCAL_NARROW_HORZ) || trig_bools.at(TriggerEnum::BitCodes::HCAL_WIDE_HORZ))
0447       {
0448         hit_threshold = 1000;
0449         waveform_hit_threshold = 1500;
0450       }
0451     }
0452   }
0453 
0454   for (int packet = packetlow; packet <= packethigh; packet++)
0455   {
0456     Packet* p = e->getPacket(packet);
0457     int zero[1] = {0};
0458     int one[1] = {1};
0459     int packet_bin = packet - packetlow + 1;
0460     if (p)
0461     {
0462       rm_packet_number[packet - packetlow]->Add(one);
0463       int packet_length[1] = {p->getLength()};
0464       rm_packet_length[packet - packetlow]->Add(packet_length);
0465 
0466       h1_packet_length->SetBinContent(packet_bin, rm_packet_length[packet - packetlow]->getMean(0));
0467 
0468       h1_packet_event->SetBinContent(packet - packetlow + 1, p->lValue(0, "CLOCK"));
0469       if (have_gl1)
0470       {
0471         long long int p_clock = p->lValue(0, "CLOCK");
0472         long long int diff = (p_clock - gl1_clock) % 65536;
0473         h_caloPack_gl1_clock_diff->Fill(packet, diff);
0474       }
0475       int nChannels = p->iValue(0, "CHANNELS");
0476       if (nChannels > m_nChannels)
0477       {
0478         return -1;  // packet is corrupted, reports too many channels
0479       }
0480       else
0481       {
0482         npacket1++;
0483         rm_packet_chans[packet - packetlow]->Add(&nChannels);
0484         h1_packet_chans->SetBinContent(packet_bin, rm_packet_chans[packet - packetlow]->getMean(0));
0485       }
0486       for (int c = 0; c < nChannels; c++)
0487       {
0488         towerNumber++;
0489 
0490         // std::vector result =  getSignal(p,c); // simple peak extraction
0491         std::vector<float> result = anaWaveform(p, c);  // full waveform fitting
0492         float signal = result.at(0);
0493         float time = result.at(1);
0494         float pedestal = result.at(2);
0495         float suppressed = result.at(result.size() - 1);
0496         float prepost = p->iValue(c, "PRE") - p->iValue(c, "POST");
0497         if (signal > 15 && signal < 15000)
0498         {
0499           energy1 += signal;
0500         }
0501 
0502         // channel mapping
0503         unsigned int key = TowerInfoDefs::encode_hcal(towerNumber - 1);
0504         unsigned int phi_bin = TowerInfoDefs::getCaloTowerPhiBin(key);
0505         unsigned int eta_bin = TowerInfoDefs::getCaloTowerEtaBin(key);
0506         int sectorNumber = phi_bin / 2 + 1;
0507         int bin = h2_hcal_mean->FindBin(eta_bin + 0.5, phi_bin + 0.5);
0508         //________________________________for this part we only want to deal with the MBD>=1 trigger
0509         if (fillhist)
0510         {
0511           if (signal > hit_threshold)
0512           {
0513             rm_vector_twrTime[towerNumber - 1]->Add(&time);
0514             rm_vector_twrhit[towerNumber - 1]->Add(one);
0515           }
0516           else
0517           {
0518             rm_vector_twrhit[towerNumber - 1]->Add(zero);
0519           }
0520           h_waveform_pedestal->Fill(pedestal);
0521 
0522           if (suppressed == 1)
0523           {
0524             pr_zsFrac_etaphi->Fill(eta_bin, phi_bin, 0);
0525           }
0526           else
0527           {
0528             pr_zsFrac_etaphi->Fill(eta_bin, phi_bin, 1);
0529           }
0530 
0531           sectorAvg[sectorNumber - 1] += signal;
0532 
0533           rm_vector_twr[towerNumber - 1]->Add(&signal);
0534 
0535           h2_hcal_mean->SetBinContent(bin, h2_hcal_mean->GetBinContent(bin) + signal);
0536           h2_hcal_rm->SetBinContent(bin, rm_vector_twrhit[towerNumber - 1]->getMean(0));
0537           h2_hcal_time->SetBinContent(bin, rm_vector_twrTime[towerNumber - 1]->getMean(0));
0538 
0539           // fill tower_rm here
0540           if (evtcnt <= historyLength * historyScaleDown)
0541           {
0542             // only fill every scaledown event
0543             if (evtcnt % historyScaleDown == 0)
0544             {
0545               h_rm_tower[eta_bin][phi_bin]->SetBinContent(evtcnt / historyScaleDown, rm_vector_twrhit[towerNumber - 1]->getMean(0));
0546             }
0547           }
0548           else
0549           {
0550             // only fill every scaledown event
0551             if (evtcnt % historyScaleDown == 0)
0552             {
0553               for (int ib = 1; ib < historyLength; ib++)
0554               {
0555                 h_rm_tower[eta_bin][phi_bin]->SetBinContent(ib, h_rm_tower[eta_bin][phi_bin]->GetBinContent(ib + 1));
0556               }
0557               h_rm_tower[eta_bin][phi_bin]->SetBinContent(historyLength, rm_vector_twrhit[towerNumber - 1]->getMean(0));
0558             }
0559           }
0560         }
0561         //_______________________________________________________end of MBD trigger requirement
0562         if (suppressed == 1)
0563         {
0564           pr_zsFrac_etaphi_all->Fill(eta_bin, phi_bin, 0);
0565         }
0566         else
0567         {
0568           pr_zsFrac_etaphi_all->Fill(eta_bin, phi_bin, 1);
0569         }
0570         // record waveform
0571         for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0572         {
0573           h_waveform_twrAvg->Fill(s, p->iValue(s, c));
0574           if (signal > waveform_hit_threshold)
0575           {
0576             h2_hcal_waveform->Fill(s, (p->iValue(s, c) - pedestal));
0577           }
0578         }
0579         if (signal > waveform_hit_threshold)
0580         {
0581           h_waveform_time->Fill(time);
0582         }
0583 
0584         if (signal > hit_threshold)
0585         {
0586           h2_hcal_hits->Fill(eta_bin + 0.5, phi_bin + 0.5);
0587           for (int itrig = 0; itrig < 64; itrig++)
0588           {
0589             if (trig_bools[itrig])
0590             {
0591               h2_hcal_hits_trig[itrig]->Fill(eta_bin + 0.5, phi_bin + 0.5);
0592             }
0593           }
0594           rm_vector_twrhit_alltrig[towerNumber - 1]->Add(one);
0595         }
0596         else
0597         {
0598           rm_vector_twrhit_alltrig[towerNumber - 1]->Add(zero);
0599         }
0600         if (prepost > 0)
0601         {
0602           p2_pre_post->Fill(eta_bin, phi_bin, prepost);
0603           p2_pre_post->Fill(eta_bin, phi_bin, -prepost);
0604         }
0605         h2_hcal_rm_alltrig->SetBinContent(bin, rm_vector_twrhit_alltrig[towerNumber - 1]->getMean(0));
0606 
0607       }  // channel loop
0608 
0609     }  // if packet good
0610     else
0611     {
0612       towerNumber += 192;
0613       rm_packet_number[packet - packetlow]->Add(zero);
0614     }
0615     h1_packet_number->SetBinContent(packet_bin, rm_packet_number[packet - packetlow]->getMean(0));
0616     delete p;
0617   }  // packet loop
0618   // if packetlow == 8001, then packetlowdiff = 7001, if packetlow == 7001, then packetlowdiff = 8001
0619   int packetlowdiff = 15002 - packetlow;
0620   int packethighdiff = 15016 - packethigh;
0621 
0622   if (npacket1 == 4)
0623   {
0624     for (int i = packetlowdiff; i <= packethighdiff; i++)
0625     {
0626       Packet* p = e->getPacket(i);
0627       if (p)
0628       {
0629         int nChannels = p->iValue(0, "CHANNELS");
0630         if (nChannels > m_nChannels)
0631         {
0632           return -1;  // packet is corrupted, reports too many channels
0633         }
0634         else
0635         {
0636           npacket2++;
0637         }
0638         for (int c = 0; c < nChannels; c++)
0639         {
0640           // std::vector result =  getSignal(p,c); // simple peak extraction
0641           std::vector<float> result = anaWaveform(p, c);  // full waveform fitting
0642           float signal = result.at(0);
0643           if (signal > 15 && signal < 15000)
0644           {
0645             energy2 += signal;
0646           }
0647         }
0648       }
0649       delete p;
0650     }
0651   }
0652   if (npacket1 == 4 && npacket2 == 4)
0653   {
0654     if (packetlow == 8001)
0655     {
0656       h2_hcal_correlation->Fill(energy1, energy2);
0657     }
0658     else
0659     {
0660       h2_hcal_correlation->Fill(energy2, energy1);
0661     }
0662   }
0663   // sector loop
0664   for (int isec = 0; isec < Nsector; isec++)
0665   {
0666     sectorAvg[isec] /= 48;
0667     h_sectorAvg_total->Fill(isec + 1, sectorAvg[isec]);
0668     rm_vector_sectAvg[isec]->Add(&sectorAvg[isec]);
0669     if (evtcnt <= historyLength * historyScaleDown)
0670     {
0671       // only fill every scaledown event
0672       if (evtcnt % historyScaleDown == 0)
0673       {
0674         h_rm_sectorAvg[isec]->SetBinContent(evtcnt / historyScaleDown, rm_vector_sectAvg[isec]->getMean(0));
0675       }
0676     }
0677     else
0678     {
0679       // only fill every scaledown event
0680       if (evtcnt % historyScaleDown == 0)
0681       {
0682         for (int ib = 1; ib < historyLength; ib++)
0683         {
0684           h_rm_sectorAvg[isec]->SetBinContent(ib, h_rm_sectorAvg[isec]->GetBinContent(ib + 1));
0685         }
0686         h_rm_sectorAvg[isec]->SetBinContent(historyLength, rm_vector_sectAvg[isec]->getMean(0));
0687       }
0688     }
0689 
0690   }  // sector loop
0691 
0692   if (fillhist) h_event->Fill(0);
0693   h_waveform_twrAvg->Scale(1. / 32. / 48.);  // average tower waveform
0694 
0695   return 0;
0696 }
0697 
0698 int HcalMon::Reset()
0699 {
0700   // reset our internal counters
0701   evtcnt = 0;
0702   idummy = 0;
0703   return 0;
0704 }