Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:20:12 

0001 #include "RecoMET/METProducers/interface/HcalNoiseInfoProducer.h"
0002 
0003 //
0004 // HcalNoiseInfoProducer.cc
0005 //
0006 //   description: Implementation of the producer for the HCAL noise information
0007 //
0008 //   author: J.P. Chou, Brown
0009 //
0010 //
0011 
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0015 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0016 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0017 #include "DataFormats/HcalRecHit/interface/HBHERecHitAuxSetter.h"
0018 #include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
0019 
0020 using namespace reco;
0021 
0022 //
0023 // constructors and destructor
0024 //
0025 
0026 HcalNoiseInfoProducer::HcalNoiseInfoProducer(const edm::ParameterSet& iConfig) : algo_(iConfig) {
0027   // set the parameters
0028   fillDigis_ = iConfig.getParameter<bool>("fillDigis");
0029   fillRecHits_ = iConfig.getParameter<bool>("fillRecHits");
0030   fillCaloTowers_ = iConfig.getParameter<bool>("fillCaloTowers");
0031   fillTracks_ = iConfig.getParameter<bool>("fillTracks");
0032   fillLaserMonitor_ = iConfig.getParameter<bool>("fillLaserMonitor");
0033 
0034   maxProblemRBXs_ = iConfig.getParameter<int>("maxProblemRBXs");
0035 
0036   maxCaloTowerIEta_ = iConfig.getParameter<int>("maxCaloTowerIEta");
0037   maxTrackEta_ = iConfig.getParameter<double>("maxTrackEta");
0038   minTrackPt_ = iConfig.getParameter<double>("minTrackPt");
0039 
0040   digiCollName_ = iConfig.getParameter<std::string>("digiCollName");
0041   recHitCollName_ = iConfig.getParameter<std::string>("recHitCollName");
0042   caloTowerCollName_ = iConfig.getParameter<std::string>("caloTowerCollName");
0043   trackCollName_ = iConfig.getParameter<std::string>("trackCollName");
0044 
0045   jetCollName_ = iConfig.getParameter<std::string>("jetCollName");
0046   maxNHF_ = iConfig.getParameter<double>("maxNHF");
0047   maxjetindex_ = iConfig.getParameter<int>("maxjetindex");
0048   jet_token_ = consumes<reco::PFJetCollection>(edm::InputTag(jetCollName_));
0049 
0050   minRecHitE_ = iConfig.getParameter<double>("minRecHitE");
0051   minLowHitE_ = iConfig.getParameter<double>("minLowHitE");
0052   minHighHitE_ = iConfig.getParameter<double>("minHighHitE");
0053 
0054   minR45HitE_ = iConfig.getParameter<double>("minR45HitE");
0055 
0056   HcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("HcalAcceptSeverityLevel");
0057   HcalRecHitFlagsToBeExcluded_ = iConfig.getParameter<std::vector<int>>("HcalRecHitFlagsToBeExcluded");
0058 
0059   // Digi threshold and time slices to use for HBHE and HF calibration digis
0060   calibdigiHBHEthreshold_ = 0;
0061   calibdigiHBHEtimeslices_ = std::vector<int>();
0062   calibdigiHFthreshold_ = 0;
0063   calibdigiHFtimeslices_ = std::vector<int>();
0064 
0065   calibdigiHBHEthreshold_ = iConfig.getParameter<double>("calibdigiHBHEthreshold");
0066   calibdigiHBHEtimeslices_ = iConfig.getParameter<std::vector<int>>("calibdigiHBHEtimeslices");
0067   calibdigiHFthreshold_ = iConfig.getParameter<double>("calibdigiHFthreshold");
0068   calibdigiHFtimeslices_ = iConfig.getParameter<std::vector<int>>("calibdigiHFtimeslices");
0069 
0070   TS4TS5EnergyThreshold_ = iConfig.getParameter<double>("TS4TS5EnergyThreshold");
0071 
0072   std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double>>("TS4TS5UpperThreshold");
0073   std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double>>("TS4TS5UpperCut");
0074   std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double>>("TS4TS5LowerThreshold");
0075   std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double>>("TS4TS5LowerCut");
0076 
0077   for (int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
0078     TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
0079   sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
0080 
0081   for (int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
0082     TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
0083   sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
0084 
0085   // if digis are filled, then rechits must also be filled
0086   if (fillDigis_ && !fillRecHits_) {
0087     fillRecHits_ = true;
0088     edm::LogWarning("HCalNoiseInfoProducer") << " forcing fillRecHits to be true if fillDigis is true.\n";
0089   }
0090 
0091   // get the fiber configuration vectors
0092   laserMonCBoxList_ = iConfig.getParameter<std::vector<int>>("laserMonCBoxList");
0093   laserMonIPhiList_ = iConfig.getParameter<std::vector<int>>("laserMonIPhiList");
0094   laserMonIEtaList_ = iConfig.getParameter<std::vector<int>>("laserMonIEtaList");
0095 
0096   // check that the vectors have the same size, if not
0097   // disable the laser monitor
0098   if (!((laserMonCBoxList_.size() == laserMonIEtaList_.size()) &&
0099         (laserMonCBoxList_.size() == laserMonIPhiList_.size()))) {
0100     edm::LogWarning("MisConfiguration") << "Must provide equally sized lists for laserMonCBoxList, laserMonIEtaList, "
0101                                            "and laserMonIPhiList.  Will not fill LaserMon\n";
0102     fillLaserMonitor_ = false;
0103   }
0104 
0105   // get the integration region with defaults
0106   laserMonitorTSStart_ = iConfig.getParameter<int>("laserMonTSStart");
0107   laserMonitorTSEnd_ = iConfig.getParameter<int>("laserMonTSEnd");
0108   laserMonitorSamples_ = iConfig.getParameter<unsigned>("laserMonSamples");
0109 
0110   adc2fC = std::vector<float>{
0111       -0.5,   0.5,    1.5,    2.5,    3.5,    4.5,    5.5,    6.5,    7.5,    8.5,    9.5,    10.5,   11.5,
0112       12.5,   13.5,   15.,    17.,    19.,    21.,    23.,    25.,    27.,    29.5,   32.5,   35.5,   38.5,
0113       42.,    46.,    50.,    54.5,   59.5,   64.5,   59.5,   64.5,   69.5,   74.5,   79.5,   84.5,   89.5,
0114       94.5,   99.5,   104.5,  109.5,  114.5,  119.5,  124.5,  129.5,  137.,   147.,   157.,   167.,   177.,
0115       187.,   197.,   209.5,  224.5,  239.5,  254.5,  272.,   292.,   312.,   334.5,  359.5,  384.5,  359.5,
0116       384.5,  409.5,  434.5,  459.5,  484.5,  509.5,  534.5,  559.5,  584.5,  609.5,  634.5,  659.5,  684.5,
0117       709.5,  747.,   797.,   847.,   897.,   947.,   997.,   1047.,  1109.5, 1184.5, 1259.5, 1334.5, 1422.,
0118       1522.,  1622.,  1734.5, 1859.5, 1984.5, 1859.5, 1984.5, 2109.5, 2234.5, 2359.5, 2484.5, 2609.5, 2734.5,
0119       2859.5, 2984.5, 3109.5, 3234.5, 3359.5, 3484.5, 3609.5, 3797.,  4047.,  4297.,  4547.,  4797.,  5047.,
0120       5297.,  5609.5, 5984.5, 6359.5, 6734.5, 7172.,  7672.,  8172.,  8734.5, 9359.5, 9984.5};
0121 
0122   // adc -> fC for qie10, for laser monitor
0123   // Taken from page 3 in
0124   // https://cms-docdb.cern.ch/cgi-bin/DocDB/RetrieveFile?docid=12570&filename=QIE10_final.pdf&version=5
0125   adc2fCHF = std::vector<float>{// - - - - - - - range 0 - - - - - - - -
0126                                 //subrange0
0127                                 1.58,
0128                                 4.73,
0129                                 7.88,
0130                                 11.0,
0131                                 14.2,
0132                                 17.3,
0133                                 20.5,
0134                                 23.6,
0135                                 26.8,
0136                                 29.9,
0137                                 33.1,
0138                                 36.2,
0139                                 39.4,
0140                                 42.5,
0141                                 45.7,
0142                                 48.8,
0143                                 //subrange1
0144                                 53.6,
0145                                 60.1,
0146                                 66.6,
0147                                 73.0,
0148                                 79.5,
0149                                 86.0,
0150                                 92.5,
0151                                 98.9,
0152                                 105,
0153                                 112,
0154                                 118,
0155                                 125,
0156                                 131,
0157                                 138,
0158                                 144,
0159                                 151,
0160                                 //subrange2
0161                                 157,
0162                                 164,
0163                                 170,
0164                                 177,
0165                                 186,
0166                                 199,
0167                                 212,
0168                                 225,
0169                                 238,
0170                                 251,
0171                                 264,
0172                                 277,
0173                                 289,
0174                                 302,
0175                                 315,
0176                                 328,
0177                                 //subrange3
0178                                 341,
0179                                 354,
0180                                 367,
0181                                 380,
0182                                 393,
0183                                 406,
0184                                 418,
0185                                 431,
0186                                 444,
0187                                 464,
0188                                 490,
0189                                 516,
0190                                 542,
0191                                 568,
0192                                 594,
0193                                 620,
0194 
0195                                 // - - - - - - - range 1 - - - - - - - -
0196                                 //subrange0
0197                                 569,
0198                                 594,
0199                                 619,
0200                                 645,
0201                                 670,
0202                                 695,
0203                                 720,
0204                                 745,
0205                                 771,
0206                                 796,
0207                                 821,
0208                                 846,
0209                                 871,
0210                                 897,
0211                                 922,
0212                                 947,
0213                                 //subrange1
0214                                 960,
0215                                 1010,
0216                                 1060,
0217                                 1120,
0218                                 1170,
0219                                 1220,
0220                                 1270,
0221                                 1320,
0222                                 1370,
0223                                 1430,
0224                                 1480,
0225                                 1530,
0226                                 1580,
0227                                 1630,
0228                                 1690,
0229                                 1740,
0230                                 //subrange2
0231                                 1790,
0232                                 1840,
0233                                 1890,
0234                                 1940,
0235                                 2020,
0236                                 2120,
0237                                 2230,
0238                                 2330,
0239                                 2430,
0240                                 2540,
0241                                 2640,
0242                                 2740,
0243                                 2850,
0244                                 2950,
0245                                 3050,
0246                                 3150,
0247                                 //subrange3
0248                                 3260,
0249                                 3360,
0250                                 3460,
0251                                 3570,
0252                                 3670,
0253                                 3770,
0254                                 3880,
0255                                 3980,
0256                                 4080,
0257                                 4240,
0258                                 4450,
0259                                 4650,
0260                                 4860,
0261                                 5070,
0262                                 5280,
0263                                 5490,
0264 
0265                                 // - - - - - - - range 2 - - - - - - - -
0266                                 //subrange0
0267                                 5080,
0268                                 5280,
0269                                 5480,
0270                                 5680,
0271                                 5880,
0272                                 6080,
0273                                 6280,
0274                                 6480,
0275                                 6680,
0276                                 6890,
0277                                 7090,
0278                                 7290,
0279                                 7490,
0280                                 7690,
0281                                 7890,
0282                                 8090,
0283                                 //subrange1
0284                                 8400,
0285                                 8810,
0286                                 9220,
0287                                 9630,
0288                                 10000,
0289                                 10400,
0290                                 10900,
0291                                 11300,
0292                                 11700,
0293                                 12100,
0294                                 12500,
0295                                 12900,
0296                                 13300,
0297                                 13700,
0298                                 14100,
0299                                 14500,
0300                                 //subrange2
0301                                 15000,
0302                                 15400,
0303                                 15800,
0304                                 16200,
0305                                 16800,
0306                                 17600,
0307                                 18400,
0308                                 19300,
0309                                 20100,
0310                                 20900,
0311                                 21700,
0312                                 22500,
0313                                 23400,
0314                                 24200,
0315                                 25000,
0316                                 25800,
0317                                 //subrange3
0318                                 26600,
0319                                 27500,
0320                                 28300,
0321                                 29100,
0322                                 29900,
0323                                 30700,
0324                                 31600,
0325                                 32400,
0326                                 33200,
0327                                 34400,
0328                                 36100,
0329                                 37700,
0330                                 39400,
0331                                 41000,
0332                                 42700,
0333                                 44300,
0334 
0335                                 // - - - - - - - range 3 - - - - - - - - -
0336                                 //subrange0
0337                                 41100,
0338                                 42700,
0339                                 44300,
0340                                 45900,
0341                                 47600,
0342                                 49200,
0343                                 50800,
0344                                 52500,
0345                                 54100,
0346                                 55700,
0347                                 57400,
0348                                 59000,
0349                                 60600,
0350                                 62200,
0351                                 63900,
0352                                 65500,
0353                                 //subrange1
0354                                 68000,
0355                                 71300,
0356                                 74700,
0357                                 78000,
0358                                 81400,
0359                                 84700,
0360                                 88000,
0361                                 91400,
0362                                 94700,
0363                                 98100,
0364                                 101000,
0365                                 105000,
0366                                 108000,
0367                                 111000,
0368                                 115000,
0369                                 118000,
0370                                 //subrange2
0371                                 121000,
0372                                 125000,
0373                                 128000,
0374                                 131000,
0375                                 137000,
0376                                 145000,
0377                                 152000,
0378                                 160000,
0379                                 168000,
0380                                 176000,
0381                                 183000,
0382                                 191000,
0383                                 199000,
0384                                 206000,
0385                                 214000,
0386                                 222000,
0387                                 //subrange3
0388                                 230000,
0389                                 237000,
0390                                 245000,
0391                                 253000,
0392                                 261000,
0393                                 268000,
0394                                 276000,
0395                                 284000,
0396                                 291000,
0397                                 302000,
0398                                 316000,
0399                                 329000,
0400                                 343000,
0401                                 356000,
0402                                 370000,
0403                                 384000};
0404 
0405   hbhedigi_token_ = consumes<HBHEDigiCollection>(edm::InputTag(digiCollName_));
0406   hcalcalibdigi_token_ = consumes<HcalCalibDigiCollection>(edm::InputTag("hcalDigis"));
0407   lasermondigi_token_ = consumes<QIE10DigiCollection>(iConfig.getParameter<edm::InputTag>("lasermonDigis"));
0408   hbherechit_token_ = consumes<HBHERecHitCollection>(edm::InputTag(recHitCollName_));
0409   calotower_token_ = consumes<CaloTowerCollection>(edm::InputTag(caloTowerCollName_));
0410   track_token_ = consumes<reco::TrackCollection>(edm::InputTag(trackCollName_));
0411   hcaltopo_token_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0412   service_token_ = esConsumes<HcalDbService, HcalDbRecord>();
0413   quality_token_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0414   severitycomputer_token_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0415   calogeometry_token_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0416 
0417   // we produce a vector of HcalNoiseRBXs
0418   produces<HcalNoiseRBXCollection>();
0419   // we also produce a noise summary
0420   produces<HcalNoiseSummary>();
0421 }
0422 
0423 HcalNoiseInfoProducer::~HcalNoiseInfoProducer() {}
0424 
0425 void HcalNoiseInfoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0426   edm::ParameterSetDescription desc;
0427   // define hit energy thesholds
0428   desc.add<double>("minRecHitE", 1.5);
0429   desc.add<double>("minLowHitE", 10.0);
0430   desc.add<double>("minHighHitE", 25.0);
0431   desc.add<double>("minR45HitE", 5.0);
0432 
0433   // define energy threshold for "problematic" cuts
0434   desc.add<double>("pMinERatio", 25.0);
0435   desc.add<double>("pMinEZeros", 5.0);
0436   desc.add<double>("pMinEEMF", 10.0);
0437 
0438   // define energy threshold for loose/tight/high level cuts
0439   desc.add<double>("minERatio", 50.0);
0440   desc.add<double>("minEZeros", 10.0);
0441   desc.add<double>("minEEMF", 50.0);
0442 
0443   // define problematic RBX
0444   desc.add<double>("pMinE", 40.0);
0445   desc.add<double>("pMinRatio", 0.75);
0446   desc.add<double>("pMaxRatio", 0.85);
0447   desc.add<int>("pMinHPDHits", 10);
0448   desc.add<int>("pMinRBXHits", 20);
0449   desc.add<int>("pMinHPDNoOtherHits", 7);
0450   desc.add<int>("pMinZeros", 4);
0451   desc.add<double>("pMinLowEHitTime", -6.0);
0452   desc.add<double>("pMaxLowEHitTime", 6.0);
0453   desc.add<double>("pMinHighEHitTime", -4.0);
0454   desc.add<double>("pMaxHighEHitTime", 5.0);
0455   desc.add<double>("pMaxHPDEMF", -0.02);
0456   desc.add<double>("pMaxRBXEMF", 0.02);
0457   desc.add<int>("pMinRBXRechitR45Count", 1);
0458   desc.add<double>("pMinRBXRechitR45Fraction", 0.1);
0459   desc.add<double>("pMinRBXRechitR45EnergyFraction", 0.1);
0460 
0461   // define loose noise cuts
0462   desc.add<double>("lMinRatio", -999.0);
0463   desc.add<double>("lMaxRatio", 999.0);
0464   desc.add<int>("lMinHPDHits", 17);
0465   desc.add<int>("lMinRBXHits", 999);
0466   desc.add<int>("lMinHPDNoOtherHits", 10);
0467   desc.add<int>("lMinZeros", 10);
0468   desc.add<double>("lMinLowEHitTime", -9999.0);
0469   desc.add<double>("lMaxLowEHitTime", 9999.0);
0470   desc.add<double>("lMinHighEHitTime", -9999.0);
0471   desc.add<double>("lMaxHighEHitTime", 9999.0);
0472 
0473   // define tight noise cuts
0474   desc.add<double>("tMinRatio", -999.0);
0475   desc.add<double>("tMaxRatio", 999.0);
0476   desc.add<int>("tMinHPDHits", 16);
0477   desc.add<int>("tMinRBXHits", 50);
0478   desc.add<int>("tMinHPDNoOtherHits", 9);
0479   desc.add<int>("tMinZeros", 8);
0480   desc.add<double>("tMinLowEHitTime", -9999.0);
0481   desc.add<double>("tMaxLowEHitTime", 9999.0);
0482   desc.add<double>("tMinHighEHitTime", -7.0);
0483   desc.add<double>("tMaxHighEHitTime", 6.0);
0484 
0485   // define high level noise cuts
0486   desc.add<double>("hlMaxHPDEMF", -9999.0);
0487   desc.add<double>("hlMaxRBXEMF", 0.01);
0488 
0489   // Calibration digi noise variables (used for finding laser noise events)
0490   desc.add<double>("calibdigiHBHEthreshold", 15)
0491       ->setComment(
0492           "minimum threshold in fC of any HBHE  \
0493               calib digi to be counted in summary");
0494   desc.add<std::vector<int>>("calibdigiHBHEtimeslices",
0495                              {
0496                                  3,
0497                                  4,
0498                                  5,
0499                                  6,
0500                              })
0501       ->setComment("time slices to use when determining charge of HBHE calib digis");
0502   desc.add<double>("calibdigiHFthreshold", -999)
0503       ->setComment("minimum threshold in fC of any HF calib digi to be counted in summary");
0504   desc.add<std::vector<int>>("calibdigiHFtimeslices",
0505                              {
0506                                  0,
0507                                  1,
0508                                  2,
0509                                  3,
0510                                  4,
0511                                  5,
0512                                  6,
0513                                  7,
0514                                  8,
0515                                  9,
0516                              })
0517       ->setComment("time slices to use when determining charge of HF calib digis");
0518 
0519   // RBX-wide TS4TS5 variable
0520   desc.add<double>("TS4TS5EnergyThreshold", 50);
0521   desc.add<std::vector<double>>("TS4TS5UpperThreshold",
0522                                 {
0523                                     70,
0524                                     90,
0525                                     100,
0526                                     400,
0527                                     4000,
0528                                 });
0529   desc.add<std::vector<double>>("TS4TS5UpperCut",
0530                                 {
0531                                     1,
0532                                     0.8,
0533                                     0.75,
0534                                     0.72,
0535                                     0.72,
0536                                 });
0537   desc.add<std::vector<double>>("TS4TS5LowerThreshold",
0538                                 {
0539                                     100,
0540                                     120,
0541                                     150,
0542                                     200,
0543                                     300,
0544                                     400,
0545                                     500,
0546                                 });
0547   desc.add<std::vector<double>>("TS4TS5LowerCut",
0548                                 {
0549                                     -1,
0550                                     -0.7,
0551                                     -0.4,
0552                                     -0.2,
0553                                     -0.08,
0554                                     0,
0555                                     0.1,
0556                                 });
0557 
0558   // rechit R45 population filter variables
0559   // this comes in groups of four: (a_Count, a_Fraction, a_EnergyFraction, const)
0560   // flag as noise if (count * a_count + fraction * a_fraction + energyfraction * a_energyfraction + const) > 0
0561   desc.add<std::vector<double>>("lRBXRecHitR45Cuts",
0562                                 {
0563                                     0.0,
0564                                     1.0,
0565                                     0.0,
0566                                     -0.5,
0567                                     0.0,
0568                                     0.0,
0569                                     1.0,
0570                                     -0.5,
0571                                 })
0572       ->setComment(
0573           "first 4 entries : equivalent to 'fraction > 0.5'  \
0574                   last 4 entries : equivalent to 'energy fraction > 0.5'");
0575   desc.add<std::vector<double>>("tRBXRecHitR45Cuts",
0576                                 {
0577                                     0.0,
0578                                     1.0,
0579                                     0.0,
0580                                     -0.2,
0581                                     0.0,
0582                                     0.0,
0583                                     1.0,
0584                                     -0.2,
0585                                 })
0586       ->setComment(
0587           "first 4 entries : equivalent to 'fraction > 0.2' \
0588                   last 4 entries : equivalent to 'energy fraction > 0.2'");
0589 
0590   // define the channels used for laser monitoring
0591   // note that the order here indicates the time order
0592   // of the channels
0593   desc.add<std::vector<int>>("laserMonCBoxList",
0594                              {
0595                                  5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
0596                              })
0597       ->setComment("time ordered list of the cBox values of laser monitor channels");
0598   desc.add<std::vector<int>>("laserMonIPhiList",
0599                              {23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0})
0600       ->setComment("time ordered list of the iPhi values of laser monitor channels");
0601   desc.add<std::vector<int>>("laserMonIEtaList",
0602                              {
0603                                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0604                              })
0605       ->setComment("time ordered list of the iEta values of laser monitor channels");
0606 
0607   // boundaries for total charge integration
0608   desc.add<int>("laserMonTSStart", 0)->setComment("lower bound of laser monitor charge integration window");
0609   desc.add<int>("laserMonTSEnd", -1)
0610       ->setComment("upper bound of laser monitor charge integration window (-1 = no bound)");
0611   desc.add<unsigned>("laserMonSamples", 4)->setComment("Number of laser monitor samples to take per channel");
0612 
0613   // what to fill
0614   desc.add<bool>("fillDigis", true);
0615   desc.add<bool>("fillRecHits", true);
0616   desc.add<bool>("fillCaloTowers", true);
0617   desc.add<bool>("fillTracks", true);
0618   desc.add<bool>("fillLaserMonitor", true);
0619 
0620   // maximum number of RBXs to fill
0621   // if you want to record all RBXs above some energy threshold,
0622   // change maxProblemRBXs to 999 and pMinE (above) to the threshold you want
0623   desc.add<int>("maxProblemRBXs", 72)
0624       ->setComment(
0625           "maximum number of RBXs to fill.  if you want to record  \
0626               all RBXs above some energy threshold,change maxProblemRBXs to  \
0627               999 and pMinE (above) to the threshold you want");
0628   ;
0629 
0630   // parameters for calculating summary variables
0631   desc.add<int>("maxCaloTowerIEta", 20);
0632   desc.add<double>("maxTrackEta", 2.0);
0633   desc.add<double>("minTrackPt", 1.0);
0634   desc.add<double>("maxNHF", 0.9);
0635   desc.add<int>("maxjetindex", 0);
0636 
0637   // collection names
0638   desc.add<std::string>("digiCollName", "hcalDigis");
0639   desc.add<std::string>("recHitCollName", "hbhereco");
0640   desc.add<std::string>("caloTowerCollName", "towerMaker");
0641   desc.add<std::string>("trackCollName", "generalTracks");
0642   desc.add<std::string>("jetCollName", "ak4PFJets");
0643   desc.add<edm::InputTag>("lasermonDigis", edm::InputTag("hcalDigis", "LASERMON"));
0644 
0645   // severity level
0646   desc.add<unsigned int>("HcalAcceptSeverityLevel", 9);
0647 
0648   // which hcal calo flags to mask
0649   // (HBHEIsolatedNoise=11, HBHEFlatNoise=12, HBHESpikeNoise=13,
0650   // HBHETriangleNoise=14, HBHETS4TS5Noise=15, HBHENegativeNoise=27)
0651   desc.add<std::vector<int>>("HcalRecHitFlagsToBeExcluded",
0652                              {
0653                                  11,
0654                                  12,
0655                                  13,
0656                                  14,
0657                                  15,
0658                                  27,
0659                              })
0660       ->setComment(
0661           "which hcal calo flags to mask (HBHEIsolatedNoise=11, \
0662                   HBHEFlatNoise=12, HBHESpikeNoise=13, \
0663                   HBHETriangleNoise=14, HBHETS4TS5Noise=15, HBHENegativeNoise=27)");
0664   ;
0665 
0666   descriptions.add("hcalnoise", desc);
0667 }
0668 
0669 //
0670 // member functions
0671 //
0672 
0673 // ------------ method called to produce the data  ------------
0674 void HcalNoiseInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0675   // this is what we're going to actually write to the EDM
0676   auto result1 = std::make_unique<HcalNoiseRBXCollection>();
0677   auto result2 = std::make_unique<HcalNoiseSummary>();
0678 
0679   // define an empty HcalNoiseRBXArray that we're going to fill
0680   HcalNoiseRBXArray rbxarray;
0681   HcalNoiseSummary& summary = *result2;
0682 
0683   // Get topology class to use later
0684   edm::ESHandle<HcalTopology> topo = iSetup.getHandle(hcaltopo_token_);
0685   theHcalTopology_ = topo.product();
0686 
0687   // fill them with the various components
0688   // digi assumes that rechit information is available
0689   if (fillRecHits_)
0690     fillrechits(iEvent, iSetup, rbxarray, summary);
0691   if (fillDigis_)
0692     filldigis(iEvent, iSetup, rbxarray, summary);
0693   if (fillCaloTowers_)
0694     fillcalotwrs(iEvent, iSetup, rbxarray, summary);
0695   if (fillTracks_)
0696     filltracks(iEvent, iSetup, summary);
0697 
0698   filljetinfo(iEvent, iSetup, summary);
0699 
0700   // select those RBXs which are interesting
0701   // also look for the highest energy RBX
0702   HcalNoiseRBXArray::iterator maxit = rbxarray.begin();
0703   double maxenergy = -999;
0704   bool maxwritten = false;
0705   for (HcalNoiseRBXArray::iterator rit = rbxarray.begin(); rit != rbxarray.end(); ++rit) {
0706     HcalNoiseRBX& rbx = (*rit);
0707     CommonHcalNoiseRBXData data(rbx,
0708                                 minRecHitE_,
0709                                 minLowHitE_,
0710                                 minHighHitE_,
0711                                 TS4TS5EnergyThreshold_,
0712                                 TS4TS5UpperCut_,
0713                                 TS4TS5LowerCut_,
0714                                 minR45HitE_);
0715 
0716     // find the highest energy rbx
0717     if (data.energy() > maxenergy) {
0718       maxenergy = data.energy();
0719       maxit = rit;
0720       maxwritten = false;
0721     }
0722 
0723     // find out if the rbx is problematic/noisy/etc.
0724     bool writerbx = algo_.isProblematic(data) || !algo_.passLooseNoiseFilter(data) ||
0725                     !algo_.passTightNoiseFilter(data) || !algo_.passHighLevelNoiseFilter(data);
0726 
0727     // fill variables in the summary object not filled elsewhere
0728     fillOtherSummaryVariables(summary, data);
0729 
0730     if (writerbx) {
0731       summary.nproblemRBXs_++;
0732       if (summary.nproblemRBXs_ <= maxProblemRBXs_) {
0733         result1->push_back(rbx);
0734         if (maxit == rit)
0735           maxwritten = true;
0736       }
0737     }
0738   }  // end loop over rbxs
0739 
0740   // if we still haven't written the maximum energy rbx, write it now
0741   if (!maxwritten) {
0742     HcalNoiseRBX& rbx = (*maxit);
0743 
0744     // add the RBX to the event
0745     result1->push_back(rbx);
0746   }
0747 
0748   // put the rbxcollection and summary into the EDM
0749   iEvent.put(std::move(result1));
0750   iEvent.put(std::move(result2));
0751 
0752   return;
0753 }
0754 
0755 // ------------ here we fill specific variables in the summary object not already accounted for earlier
0756 void HcalNoiseInfoProducer::fillOtherSummaryVariables(HcalNoiseSummary& summary,
0757                                                       const CommonHcalNoiseRBXData& data) const {
0758   // charge ratio
0759   if (algo_.passRatioThreshold(data) && data.validRatio()) {
0760     if (data.ratio() < summary.minE2Over10TS()) {
0761       summary.mine2ts_ = data.e2ts();
0762       summary.mine10ts_ = data.e10ts();
0763     }
0764     if (data.ratio() > summary.maxE2Over10TS()) {
0765       summary.maxe2ts_ = data.e2ts();
0766       summary.maxe10ts_ = data.e10ts();
0767     }
0768   }
0769 
0770   // ADC zeros
0771   if (algo_.passZerosThreshold(data)) {
0772     if (data.numZeros() > summary.maxZeros()) {
0773       summary.maxzeros_ = data.numZeros();
0774     }
0775   }
0776 
0777   // hits count
0778   if (data.numHPDHits() > summary.maxHPDHits()) {
0779     summary.maxhpdhits_ = data.numHPDHits();
0780   }
0781   if (data.numRBXHits() > summary.maxRBXHits()) {
0782     summary.maxrbxhits_ = data.numRBXHits();
0783   }
0784   if (data.numHPDNoOtherHits() > summary.maxHPDNoOtherHits()) {
0785     summary.maxhpdhitsnoother_ = data.numHPDNoOtherHits();
0786   }
0787 
0788   // TS4TS5
0789   if (data.PassTS4TS5() == false)
0790     summary.hasBadRBXTS4TS5_ = true;
0791 
0792   if (algo_.passLooseRBXRechitR45(data) == false)
0793     summary.hasBadRBXRechitR45Loose_ = true;
0794   if (algo_.passTightRBXRechitR45(data) == false)
0795     summary.hasBadRBXRechitR45Tight_ = true;
0796 
0797   // hit timing
0798   if (data.minLowEHitTime() < summary.min10GeVHitTime()) {
0799     summary.min10_ = data.minLowEHitTime();
0800   }
0801   if (data.maxLowEHitTime() > summary.max10GeVHitTime()) {
0802     summary.max10_ = data.maxLowEHitTime();
0803   }
0804   summary.rms10_ += data.lowEHitTimeSqrd();
0805   summary.cnthit10_ += data.numLowEHits();
0806   if (data.minHighEHitTime() < summary.min25GeVHitTime()) {
0807     summary.min25_ = data.minHighEHitTime();
0808   }
0809   if (data.maxHighEHitTime() > summary.max25GeVHitTime()) {
0810     summary.max25_ = data.maxHighEHitTime();
0811   }
0812   summary.rms25_ += data.highEHitTimeSqrd();
0813   summary.cnthit25_ += data.numHighEHits();
0814 
0815   // EMF
0816   if (algo_.passEMFThreshold(data)) {
0817     if (summary.minHPDEMF() > data.HPDEMF()) {
0818       summary.minhpdemf_ = data.HPDEMF();
0819     }
0820     if (summary.minRBXEMF() > data.RBXEMF()) {
0821       summary.minrbxemf_ = data.RBXEMF();
0822     }
0823   }
0824 
0825   // summary flag
0826   if (!algo_.passLooseRatio(data))
0827     summary.filterstatus_ |= 0x1;
0828   if (!algo_.passLooseHits(data))
0829     summary.filterstatus_ |= 0x2;
0830   if (!algo_.passLooseZeros(data))
0831     summary.filterstatus_ |= 0x4;
0832   if (!algo_.passLooseTiming(data))
0833     summary.filterstatus_ |= 0x8;
0834 
0835   if (!algo_.passTightRatio(data))
0836     summary.filterstatus_ |= 0x100;
0837   if (!algo_.passTightHits(data))
0838     summary.filterstatus_ |= 0x200;
0839   if (!algo_.passTightZeros(data))
0840     summary.filterstatus_ |= 0x400;
0841   if (!algo_.passTightTiming(data))
0842     summary.filterstatus_ |= 0x800;
0843 
0844   if (!algo_.passHighLevelNoiseFilter(data))
0845     summary.filterstatus_ |= 0x10000;
0846 
0847   // summary refvectors
0848   JoinCaloTowerRefVectorsWithoutDuplicates join;
0849   if (!algo_.passLooseNoiseFilter(data))
0850     join(summary.loosenoisetwrs_, data.rbxTowers());
0851   if (!algo_.passTightNoiseFilter(data))
0852     join(summary.tightnoisetwrs_, data.rbxTowers());
0853   if (!algo_.passHighLevelNoiseFilter(data))
0854     join(summary.hlnoisetwrs_, data.rbxTowers());
0855 
0856   return;
0857 }
0858 
0859 // ------------ fill the array with digi information
0860 void HcalNoiseInfoProducer::filldigis(edm::Event& iEvent,
0861                                       const edm::EventSetup& iSetup,
0862                                       HcalNoiseRBXArray& array,
0863                                       HcalNoiseSummary& summary) {
0864   // Some initialization
0865   totalCalibCharge = 0;
0866   totalLasmonCharge = 0;
0867 
0868   // Starting with this version (updated by Jeff Temple on Dec. 6, 2012), the "TS45" names in the variables are mis-nomers.  The actual time slices used are determined from the digiTimeSlices_ variable, which may not be limited to only time slices 4 and 5.  For now, "TS45" name kept, because that is what is used in HcalNoiseSummary object (in GetCalibCountTS45, etc.).  Likewise, the charge value in 'gt15' is now configurable, though the name remains the same.  For HBHE, we track both the number of calibration channels (NcalibTS45) and the number of calibration channels above threshold (NcalibTS45gt15).  For HF, we track only the number of channels above the given threshold in the given time window (NcalibHFgtX).  Default for HF in 2012 is to use the full time sample with effectively no threshold (threshold=-999)
0869   int NcalibTS45 = 0;
0870   int NcalibTS45gt15 = 0;
0871   int NcalibHFgtX = 0;
0872 
0873   double chargecalibTS45 = 0;
0874   double chargecalibgt15TS45 = 0;
0875 
0876   // get the conditions and channel quality
0877   edm::ESHandle<HcalDbService> conditions = iSetup.getHandle(service_token_);
0878   edm::ESHandle<HcalChannelQuality> qualhandle = iSetup.getHandle(quality_token_);
0879   const HcalChannelQuality* myqual = qualhandle.product();
0880   edm::ESHandle<HcalSeverityLevelComputer> mycomputer = iSetup.getHandle(severitycomputer_token_);
0881   const HcalSeverityLevelComputer* mySeverity = mycomputer.product();
0882 
0883   // get the digis
0884   edm::Handle<HBHEDigiCollection> handle;
0885   //  iEvent.getByLabel(digiCollName_, handle);
0886   iEvent.getByToken(hbhedigi_token_, handle);
0887 
0888   if (!handle.isValid()) {
0889     throw edm::Exception(edm::errors::ProductNotFound)
0890         << " could not find HBHEDigiCollection named " << digiCollName_ << "\n.";
0891     return;
0892   }
0893 
0894   // loop over all of the digi information
0895   for (HBHEDigiCollection::const_iterator it = handle->begin(); it != handle->end(); ++it) {
0896     const HBHEDataFrame& digi = (*it);
0897     HcalDetId cell = digi.id();
0898     DetId detcell = (DetId)cell;
0899 
0900     // check on cells to be ignored and dropped
0901     const HcalChannelStatus* mydigistatus = myqual->getValues(detcell.rawId());
0902     if (mySeverity->dropChannel(mydigistatus->getValue()))
0903       continue;
0904     if (digi.zsMarkAndPass())
0905       continue;
0906     // Drop if exclude bit set
0907     if (mydigistatus->isBitSet(HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummary))
0908       continue;
0909 
0910     // get the calibrations and coder
0911     const HcalCalibrations& calibrations = conditions->getHcalCalibrations(cell);
0912     const HcalQIECoder* channelCoder = conditions->getHcalCoder(cell);
0913     const HcalQIEShape* shape = conditions->getHcalShape(channelCoder);
0914     HcalCoderDb coder(*channelCoder, *shape);
0915 
0916     // match the digi to an rbx and hpd
0917     HcalNoiseRBX& rbx = (*array.findRBX(digi));
0918     HcalNoiseHPD& hpd = (*array.findHPD(digi));
0919 
0920     // determine if the digi is one the highest energy hits in the HPD
0921     // this works because the rechits are sorted by energy (see fillrechits() below)
0922     bool isBig = false, isBig5 = false, isRBX = false;
0923     int counter = 0;
0924     edm::RefVector<HBHERecHitCollection>& rechits = hpd.rechits_;
0925     for (edm::RefVector<HBHERecHitCollection>::const_iterator rit = rechits.begin(); rit != rechits.end();
0926          ++rit, ++counter) {
0927       const HcalDetId& detid = (*rit)->idFront();
0928       if (DetId(detid) == digi.id()) {
0929         if (counter == 0)
0930           isBig = isBig5 = true;  // digi is also the highest energy rechit
0931         if (counter < 5)
0932           isBig5 = true;  // digi is one of 5 highest energy rechits
0933         isRBX = true;
0934       }
0935     }
0936 
0937     // loop over each of the digi's time slices
0938     int totalzeros = 0;
0939     CaloSamples tool;
0940     coder.adc2fC(digi, tool);
0941     for (int ts = 0; ts < tool.size(); ++ts) {
0942       // count zero's
0943       if (digi[ts].adc() == 0) {
0944         ++hpd.totalZeros_;
0945         ++totalzeros;
0946       }
0947 
0948       // get the fC's
0949       double corrfc = tool[ts] - calibrations.pedestal(digi[ts].capid());
0950 
0951       // fill the relevant digi arrays
0952       if (isBig)
0953         hpd.bigCharge_[ts] += corrfc;
0954       if (isBig5)
0955         hpd.big5Charge_[ts] += corrfc;
0956       if (isRBX)
0957         rbx.allCharge_[ts] += corrfc;
0958     }
0959 
0960     // record the maximum number of zero's found
0961     if (totalzeros > hpd.maxZeros_)
0962       hpd.maxZeros_ = totalzeros;
0963   }
0964 
0965   // get the calibration digis
0966   edm::Handle<HcalCalibDigiCollection> hCalib;
0967   //  iEvent.getByLabel("hcalDigis", hCalib);
0968   iEvent.getByToken(hcalcalibdigi_token_, hCalib);
0969 
0970   // get the lasermon digis
0971   edm::Handle<QIE10DigiCollection> hLasermon;
0972   iEvent.getByToken(lasermondigi_token_, hLasermon);
0973 
0974   // get total charge in calibration channels
0975   if (hCalib.isValid() == true) {
0976     for (HcalCalibDigiCollection::const_iterator digi = hCalib->begin(); digi != hCalib->end(); digi++) {
0977       if (digi->id().hcalSubdet() == 0)
0978         continue;
0979 
0980       for (unsigned i = 0; i < (unsigned)digi->size(); i++)
0981         totalCalibCharge = totalCalibCharge + adc2fC[digi->sample(i).adc() & 0xff];
0982 
0983       HcalCalibDetId myid = (HcalCalibDetId)digi->id();
0984       if (myid.calibFlavor() == HcalCalibDetId::HOCrosstalk)
0985         continue;  // ignore HOCrosstalk channels
0986       if (digi->zsMarkAndPass())
0987         continue;  // skip "mark-and-pass" channels when computing charge in calib channels
0988 
0989       if (digi->id().hcalSubdet() == HcalForward)  // check HF
0990       {
0991         double sumChargeHF = 0;
0992         for (unsigned int i = 0; i < calibdigiHFtimeslices_.size(); ++i) {
0993           // skip unphysical time slices
0994           if (calibdigiHFtimeslices_[i] < 0 || calibdigiHFtimeslices_[i] > digi->size())
0995             continue;
0996           sumChargeHF += adc2fC[digi->sample(calibdigiHFtimeslices_[i]).adc() & 0xff];
0997         }
0998         if (sumChargeHF > calibdigiHFthreshold_)
0999           ++NcalibHFgtX;
1000       }                                                                                         // end of HF check
1001       else if (digi->id().hcalSubdet() == HcalBarrel || digi->id().hcalSubdet() == HcalEndcap)  // now check HBHE
1002       {
1003         double sumChargeHBHE = 0;
1004         for (unsigned int i = 0; i < calibdigiHBHEtimeslices_.size(); ++i) {
1005           // skip unphysical time slices
1006           if (calibdigiHBHEtimeslices_[i] < 0 || calibdigiHBHEtimeslices_[i] > digi->size())
1007             continue;
1008           sumChargeHBHE += adc2fC[digi->sample(calibdigiHBHEtimeslices_[i]).adc() & 0xff];
1009         }
1010         ++NcalibTS45;
1011         chargecalibTS45 += sumChargeHBHE;
1012         if (sumChargeHBHE > calibdigiHBHEthreshold_) {
1013           ++NcalibTS45gt15;
1014           chargecalibgt15TS45 += sumChargeHBHE;
1015         }
1016       }  // end of HBHE check
1017     }    // loop on HcalCalibDigiCollection
1018 
1019   }  // if (hCalib.isValid()==true)
1020   if (fillLaserMonitor_ && (hLasermon.isValid() == true)) {
1021     int icombts = -1;
1022     float max_charge = 0;
1023     int max_ts = -1;
1024     std::vector<float> comb_charge;
1025 
1026     unsigned nch = laserMonCBoxList_.size();
1027     // loop over channels in the order provided
1028     for (unsigned ich = 0; ich < nch; ++ich) {
1029       int cboxch = laserMonCBoxList_[ich];
1030       int iphi = laserMonIPhiList_[ich];
1031       int ieta = laserMonIEtaList_[ich];
1032 
1033       // loop over digis, find the digi that matches this channel
1034       for (const QIE10DataFrame df : (*hLasermon)) {
1035         HcalCalibDetId calibId(df.id());
1036 
1037         int ch_cboxch = calibId.cboxChannel();
1038         int ch_iphi = calibId.iphi();
1039         int ch_ieta = calibId.ieta();
1040 
1041         if (cboxch == ch_cboxch && iphi == ch_iphi && ieta == ch_ieta) {
1042           unsigned ts_size = df.samples();
1043 
1044           // loop over time slices
1045           for (unsigned its = 0; its < ts_size; ++its) {
1046             // if we are on the last channel, use all the data
1047             // otherwise only take the unique samples
1048             if (((ich + 1) < nch) && its >= laserMonitorSamples_)
1049               continue;
1050 
1051             bool ok = df[its].ok();
1052             int adc = df[its].adc();
1053 
1054             icombts++;
1055             // apply integration limits
1056             if (icombts < laserMonitorTSStart_)
1057               continue;
1058             if (laserMonitorTSEnd_ > 0 && icombts > laserMonitorTSEnd_)
1059               continue;
1060 
1061             if (ok && adc >= 0) {  // protection against QIE reset or bad ADC values
1062 
1063               float charge = adc2fCHF[adc];
1064               if (charge > max_charge) {
1065                 max_charge = charge;
1066                 max_ts = icombts;
1067               }
1068 
1069               comb_charge.push_back(charge);
1070             }  // if( ok && adc >= 0 )
1071           }    // loop over time slices
1072         }      // if( cboxch == ch_cboxch && iphi == ch_iphi && ieta == ch_ieta )
1073       }        // loop over digi collection
1074     }          // loop over channel list
1075 
1076     // do not continue with the calculation
1077     // if the vector was not filled
1078     if (comb_charge.empty()) {
1079       totalLasmonCharge = -1;
1080     } else {
1081       // integrate from +- 3 TS around the time sample
1082       // having the maximum charge
1083       int start_ts = max_ts - 3;
1084       int end_ts = max_ts + 3;
1085 
1086       // Change the integration limits
1087       // if outside of the range
1088       if (start_ts < 0)
1089         start_ts = 0;
1090       if (end_ts >= int(comb_charge.size()))
1091         end_ts = comb_charge.size() - 1;
1092 
1093       for (int isum = start_ts; isum <= end_ts; ++isum) {
1094         totalLasmonCharge += comb_charge[isum];
1095       }
1096     }
1097   }  // if( fillLaserMonitor_  && (hLasermon.isValid() == true) )
1098 
1099   summary.calibCharge_ = totalCalibCharge;
1100   summary.lasmonCharge_ = totalLasmonCharge;
1101   summary.calibCountTS45_ = NcalibTS45;
1102   summary.calibCountgt15TS45_ = NcalibTS45gt15;
1103   summary.calibChargeTS45_ = chargecalibTS45;
1104   summary.calibChargegt15TS45_ = chargecalibgt15TS45;
1105   summary.calibCountHF_ = NcalibHFgtX;
1106 
1107   return;
1108 }
1109 
1110 // ------------ fill the array with rec hit information
1111 void HcalNoiseInfoProducer::fillrechits(edm::Event& iEvent,
1112                                         const edm::EventSetup& iSetup,
1113                                         HcalNoiseRBXArray& array,
1114                                         HcalNoiseSummary& summary) const {
1115   // get the HCAL channel status map
1116   edm::ESHandle<HcalChannelQuality> hcalChStatus = iSetup.getHandle(quality_token_);
1117   const HcalChannelQuality* dbHcalChStatus = hcalChStatus.product();
1118 
1119   // get the severity level computer
1120   edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl = iSetup.getHandle(severitycomputer_token_);
1121   const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
1122 
1123   // get the calo geometry
1124   edm::ESHandle<CaloGeometry> pG = iSetup.getHandle(calogeometry_token_);
1125   const CaloGeometry* geo = pG.product();
1126 
1127   // get the rechits
1128   edm::Handle<HBHERecHitCollection> handle;
1129   //  iEvent.getByLabel(recHitCollName_, handle);
1130   iEvent.getByToken(hbherechit_token_, handle);
1131 
1132   if (!handle.isValid()) {
1133     throw edm::Exception(edm::errors::ProductNotFound)
1134         << " could not find HBHERecHitCollection named " << recHitCollName_ << "\n.";
1135     return;
1136   }
1137 
1138   summary.rechitCount_ = 0;
1139   summary.rechitCount15_ = 0;
1140   summary.rechitEnergy_ = 0;
1141   summary.rechitEnergy15_ = 0;
1142 
1143   summary.hitsInLaserRegion_ = 0;
1144   summary.hitsInNonLaserRegion_ = 0;
1145   summary.energyInLaserRegion_ = 0;
1146   summary.energyInNonLaserRegion_ = 0;
1147 
1148   // loop over all of the rechit information
1149   for (HBHERecHitCollection::const_iterator it = handle->begin(); it != handle->end(); ++it) {
1150     const HBHERecHit& rechit = (*it);
1151 
1152     // skip bad rechits (other than those flagged by the isolated noise, triangle, flat, and spike algorithms)
1153     const DetId id = rechit.idFront();
1154 
1155     uint32_t recHitFlag = rechit.flags();
1156     uint32_t isolbitset = (1 << HcalCaloFlagLabels::HBHEIsolatedNoise);
1157     uint32_t flatbitset = (1 << HcalCaloFlagLabels::HBHEFlatNoise);
1158     uint32_t spikebitset = (1 << HcalCaloFlagLabels::HBHESpikeNoise);
1159     uint32_t trianglebitset = (1 << HcalCaloFlagLabels::HBHETriangleNoise);
1160     uint32_t ts4ts5bitset = (1 << HcalCaloFlagLabels::HBHETS4TS5Noise);
1161     uint32_t negativebitset = (1 << HcalCaloFlagLabels::HBHENegativeNoise);
1162     for (unsigned int i = 0; i < HcalRecHitFlagsToBeExcluded_.size(); i++) {
1163       uint32_t bitset = (1 << HcalRecHitFlagsToBeExcluded_[i]);
1164       recHitFlag = (recHitFlag & bitset) ? recHitFlag - bitset : recHitFlag;
1165     }
1166     const HcalChannelStatus* dbStatus = dbHcalChStatus->getValues(id);
1167 
1168     // Ignore rechit if exclude bit set, regardless of severity of other bits
1169     if (dbStatus->isBitSet(HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummary))
1170       continue;
1171 
1172     int severityLevel = hcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatus->getValue());
1173     bool isRecovered = hcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
1174     if (severityLevel != 0 && !isRecovered && severityLevel > static_cast<int>(HcalAcceptSeverityLevel_))
1175       continue;
1176 
1177     // do some rechit counting and energies
1178     summary.rechitCount_ = summary.rechitCount_ + 1;
1179     summary.rechitEnergy_ = summary.rechitEnergy_ + rechit.eraw();
1180     if (dbStatus->isBitSet(
1181             HcalChannelStatus::
1182                 HcalBadLaserSignal))  // hit comes from a region where no laser calibration pulse is normally seen
1183     {
1184       ++summary.hitsInNonLaserRegion_;
1185       summary.energyInNonLaserRegion_ += rechit.eraw();
1186     } else  // hit comes from region where laser calibration pulse is seen
1187     {
1188       ++summary.hitsInLaserRegion_;
1189       summary.energyInLaserRegion_ += rechit.eraw();
1190     }
1191 
1192     if (rechit.eraw() > 1.5) {
1193       summary.rechitCount15_ = summary.rechitCount15_ + 1;
1194       summary.rechitEnergy15_ = summary.rechitEnergy15_ + rechit.eraw();
1195     }
1196 
1197     // Exclude uncollapsed QIE11 channels
1198     if (CaloRecHitAuxSetter::getBit(rechit.auxPhase1(), HBHERecHitAuxSetter::OFF_TDC_TIME) &&
1199         !CaloRecHitAuxSetter::getBit(rechit.auxPhase1(), HBHERecHitAuxSetter::OFF_COMBINED))
1200       continue;
1201 
1202     // if it was ID'd as isolated noise, update the summary object
1203     if (rechit.flags() & isolbitset) {
1204       summary.nisolnoise_++;
1205       summary.isolnoisee_ += rechit.eraw();
1206       GlobalPoint gp = geo->getPosition(rechit.id());
1207       double et = rechit.eraw() * gp.perp() / gp.mag();
1208       summary.isolnoiseet_ += et;
1209     }
1210 
1211     if (rechit.flags() & flatbitset) {
1212       summary.nflatnoise_++;
1213       summary.flatnoisee_ += rechit.eraw();
1214       GlobalPoint gp = geo->getPosition(rechit.id());
1215       double et = rechit.eraw() * gp.perp() / gp.mag();
1216       summary.flatnoiseet_ += et;
1217     }
1218 
1219     if (rechit.flags() & spikebitset) {
1220       summary.nspikenoise_++;
1221       summary.spikenoisee_ += rechit.eraw();
1222       GlobalPoint gp = geo->getPosition(rechit.id());
1223       double et = rechit.eraw() * gp.perp() / gp.mag();
1224       summary.spikenoiseet_ += et;
1225     }
1226 
1227     if (rechit.flags() & trianglebitset) {
1228       summary.ntrianglenoise_++;
1229       summary.trianglenoisee_ += rechit.eraw();
1230       GlobalPoint gp = geo->getPosition(rechit.id());
1231       double et = rechit.eraw() * gp.perp() / gp.mag();
1232       summary.trianglenoiseet_ += et;
1233     }
1234 
1235     if (rechit.flags() & ts4ts5bitset) {
1236       // only add to TS4TS5 if the bit is not marked as "HcalCellExcludeFromHBHENoiseSummaryR45"
1237       if (not dbStatus->isBitSet(HcalChannelStatus::HcalCellExcludeFromHBHENoiseSummaryR45)) {
1238         summary.nts4ts5noise_++;
1239         summary.ts4ts5noisee_ += rechit.eraw();
1240         GlobalPoint gp = geo->getPosition(rechit.id());
1241         double et = rechit.eraw() * gp.perp() / gp.mag();
1242         summary.ts4ts5noiseet_ += et;
1243       }
1244     }
1245 
1246     if (rechit.flags() & negativebitset) {
1247       summary.nnegativenoise_++;
1248       summary.negativenoisee_ += rechit.eraw();
1249       GlobalPoint gp = geo->getPosition(rechit.id());
1250       double et = rechit.eraw() * gp.perp() / gp.mag();
1251       summary.negativenoiseet_ += et;
1252     }
1253 
1254     // find the hpd that the rechit is in
1255     HcalNoiseHPD& hpd = (*array.findHPD(rechit));
1256 
1257     // create a persistent reference to the rechit
1258     edm::Ref<HBHERecHitCollection> myRef(handle, it - handle->begin());
1259 
1260     // store it in a place so that it remains sorted by energy
1261     hpd.refrechitset_.insert(myRef);
1262 
1263   }  // end loop over rechits
1264 
1265   // loop over all HPDs and transfer the information from refrechitset_ to rechits_;
1266   for (HcalNoiseRBXArray::iterator rbxit = array.begin(); rbxit != array.end(); ++rbxit) {
1267     for (std::vector<HcalNoiseHPD>::iterator hpdit = rbxit->hpds_.begin(); hpdit != rbxit->hpds_.end(); ++hpdit) {
1268       // loop over all of the entries in the set and add them to rechits_
1269       for (std::set<edm::Ref<HBHERecHitCollection>, RefHBHERecHitEnergyComparison>::const_iterator it =
1270                hpdit->refrechitset_.begin();
1271            it != hpdit->refrechitset_.end();
1272            ++it) {
1273         hpdit->rechits_.push_back(*it);
1274       }
1275     }
1276   }
1277   // now the rechits in all the HPDs are sorted by energy!
1278 
1279   return;
1280 }
1281 
1282 // ------------ fill the array with calo tower information
1283 void HcalNoiseInfoProducer::fillcalotwrs(edm::Event& iEvent,
1284                                          const edm::EventSetup& iSetup,
1285                                          HcalNoiseRBXArray& array,
1286                                          HcalNoiseSummary& summary) const {
1287   // get the calotowers
1288   edm::Handle<CaloTowerCollection> handle;
1289   //  iEvent.getByLabel(caloTowerCollName_, handle);
1290   iEvent.getByToken(calotower_token_, handle);
1291 
1292   if (!handle.isValid()) {
1293     throw edm::Exception(edm::errors::ProductNotFound)
1294         << " could not find CaloTowerCollection named " << caloTowerCollName_ << "\n.";
1295     return;
1296   }
1297 
1298   summary.emenergy_ = summary.hadenergy_ = 0.0;
1299 
1300   // loop over all of the calotower information
1301   for (CaloTowerCollection::const_iterator it = handle->begin(); it != handle->end(); ++it) {
1302     const CaloTower& twr = (*it);
1303 
1304     // create a persistent reference to the tower
1305     edm::Ref<CaloTowerCollection> myRef(handle, it - handle->begin());
1306 
1307     // get all of the hpd's that are pointed to by the calotower
1308     std::vector<std::vector<HcalNoiseHPD>::iterator> hpditervec;
1309     array.findHPD(twr, hpditervec);
1310 
1311     // loop over the hpd's and add the reference to the RefVectors
1312     for (std::vector<std::vector<HcalNoiseHPD>::iterator>::iterator it = hpditervec.begin(); it != hpditervec.end();
1313          ++it)
1314       (*it)->calotowers_.push_back(myRef);
1315 
1316     // skip over anything with |ieta|>maxCaloTowerIEta
1317     if (twr.ietaAbs() > maxCaloTowerIEta_) {
1318       summary.emenergy_ += twr.emEnergy();
1319       summary.hadenergy_ += twr.hadEnergy();
1320     }
1321   }
1322 
1323   return;
1324 }
1325 
1326 // ------------ fill the summary info from jets
1327 void HcalNoiseInfoProducer::filljetinfo(edm::Event& iEvent,
1328                                         const edm::EventSetup& iSetup,
1329                                         HcalNoiseSummary& summary) const {
1330   bool goodJetFoundInLowBVRegion = false;  // checks whether a jet is in
1331                                            // a low BV region, where false
1332                                            // noise flagging rate is higher.
1333   if (!jetCollName_.empty()) {
1334     edm::Handle<reco::PFJetCollection> pfjet_h;
1335     iEvent.getByToken(jet_token_, pfjet_h);
1336 
1337     if (pfjet_h.isValid()) {
1338       int jetindex = 0;
1339       for (reco::PFJetCollection::const_iterator jet = pfjet_h->begin(); jet != pfjet_h->end(); ++jet) {
1340         if (jetindex > maxjetindex_)
1341           break;  // only look at jets with
1342                   // indices up to maxjetindex_
1343 
1344         // Check whether jet is in low-BV region (0<eta<1.4, -1.8<phi<-1.4)
1345         if (jet->eta() > 0.0 && jet->eta() < 1.4 && jet->phi() > -1.8 && jet->phi() < -1.4) {
1346           // Look for a good jet in low BV region;
1347           // if found, we will keep event
1348           if (maxNHF_ < 0.0 || jet->neutralHadronEnergyFraction() < maxNHF_) {
1349             goodJetFoundInLowBVRegion = true;
1350             break;
1351           }
1352         }
1353         ++jetindex;
1354       }
1355     }
1356   }
1357 
1358   summary.goodJetFoundInLowBVRegion_ = goodJetFoundInLowBVRegion;
1359 }
1360 
1361 // ------------ fill the summary with track information
1362 void HcalNoiseInfoProducer::filltracks(edm::Event& iEvent,
1363                                        const edm::EventSetup& iSetup,
1364                                        HcalNoiseSummary& summary) const {
1365   edm::Handle<reco::TrackCollection> handle;
1366   //  iEvent.getByLabel(trackCollName_, handle);
1367   iEvent.getByToken(track_token_, handle);
1368 
1369   // don't throw exception, just return quietly
1370   if (!handle.isValid()) {
1371     //    throw edm::Exception(edm::errors::ProductNotFound)
1372     //      << " could not find trackCollection named " << trackCollName_ << "\n.";
1373     return;
1374   }
1375 
1376   summary.trackenergy_ = 0.0;
1377   for (reco::TrackCollection::const_iterator iTrack = handle->begin(); iTrack != handle->end(); ++iTrack) {
1378     const reco::Track& trk = *iTrack;
1379     if (trk.pt() < minTrackPt_ || fabs(trk.eta()) > maxTrackEta_)
1380       continue;
1381 
1382     summary.trackenergy_ += trk.p();
1383   }
1384 
1385   return;
1386 }
1387 
1388 //define this as a plug-in
1389 DEFINE_FWK_MODULE(HcalNoiseInfoProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!