METReco/interface/AnomalousECALVariables.h

0001 #ifndef ANOMALOUSECALVARIABLES_H_
0002 #define ANOMALOUSECALVARIABLES_H_
0003 //DataFormats/AnomalousEcalDataFormats/interface/AnomalousECALVariables.h
0004 // system include files
0005 #include <memory>
0006
0007 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0008 #include "DataFormats/METReco/interface/BoundaryInformation.h"
0009
0010 //using namespace edm;
0011 //using namespace std;
0012 /*
0013  * This class summarizes the information about the boundary energy calculated in EcalAnomalousEventFilter:
0014  * 1. next to ECAL border/gap
0015  * 2. next to masked ECAL channels: for each dead area with boundary energy above a threshold defined in the filter
0016  * the vector 'v_enDeadNeighbours_EB' or 'v_enDeadNeighbours_EE' is filled with the calculated boundary energy.
0017  * The determined size of the corresponding cluster is filled in v_enDeadNeighboursNoCells_EB/EE accordingly.
0018  *
0019  */
0020 class AnomalousECALVariables {
0021 public:
0022   AnomalousECALVariables() {
0023     //energy next to ECAL Gap
0024     v_enNeighboursGap_EB.reserve(50);
0025     v_enNeighboursGap_EE.reserve(50);
0026     v_enNeighboursGap_EB.clear();
0027     v_enNeighboursGap_EE.clear();
0028
0029     //energy around dead cells
0030     v_boundaryInfoDeadCells_EB = std::vector<BoundaryInformation>();
0031     v_boundaryInfoDeadCells_EE = std::vector<BoundaryInformation>();
0032     v_boundaryInfoDeadCells_EB.reserve(50);
0033     v_boundaryInfoDeadCells_EE.reserve(50);
0034     v_boundaryInfoDeadCells_EB.clear();
0035     v_boundaryInfoDeadCells_EE.clear();
0036   };
0037
0038   AnomalousECALVariables(const std::vector<BoundaryInformation>& p_enNeighboursGap_EB,
0039                          const std::vector<BoundaryInformation>& p_enNeighboursGap_EE,
0040                          const std::vector<BoundaryInformation>& p_boundaryInfoDeadCells_EB,
0041                          const std::vector<BoundaryInformation>& p_boundaryInfoDeadCells_EE) {
0042     v_boundaryInfoDeadCells_EB = std::vector<BoundaryInformation>();
0043     v_boundaryInfoDeadCells_EE = std::vector<BoundaryInformation>();
0044     v_boundaryInfoDeadCells_EB.reserve(50);
0045     v_boundaryInfoDeadCells_EE.reserve(50);
0046     v_boundaryInfoDeadCells_EB.clear();
0047     v_boundaryInfoDeadCells_EE.clear();
0048     v_boundaryInfoDeadCells_EB = p_boundaryInfoDeadCells_EB;
0049     v_boundaryInfoDeadCells_EE = p_boundaryInfoDeadCells_EE;
0050
0051     v_enNeighboursGap_EB = p_enNeighboursGap_EB;
0052     v_enNeighboursGap_EE = p_enNeighboursGap_EE;
0053   };
0054
0055   ~AnomalousECALVariables() {
0056     //cout<<"destructor AnomalousECAL"<<endl;
0057     v_enNeighboursGap_EB.clear();
0058     v_enNeighboursGap_EE.clear();
0059     v_boundaryInfoDeadCells_EB.clear();
0060     v_boundaryInfoDeadCells_EE.clear();
0061   };
0062
0063   //returns true if at least 1 dead cell area was found in EcalAnomalousEventFilter with
0064   //boundary energy above threshold
0065   //Note: no sense to change this cut BELOW the threshold given in EcalAnomalousEventFilter
0066
0067   bool isDeadEcalCluster(double maxBoundaryEnergy = 10,
0068                          const std::vector<int>& limitDeadCellToChannelStatusEB = std::vector<int>(),
0069                          const std::vector<int>& limitDeadCellToChannelStatusEE = std::vector<int>()) const {
0070     float highestEnergyDepositAroundDeadCell = 0;
0071
0072     for (int i = 0; i < (int)v_boundaryInfoDeadCells_EB.size(); ++i) {
0073       BoundaryInformation bInfo = v_boundaryInfoDeadCells_EB[i];
0074
0075       //check if channel limitation rejectsbInfo
0076       bool passChannelLimitation = false;
0077       std::vector<int> status = bInfo.channelStatus;
0078
0079       for (int cs = 0; cs < (int)limitDeadCellToChannelStatusEB.size(); ++cs) {
0080         int channelAllowed = limitDeadCellToChannelStatusEB[cs];
0081
0082         for (std::vector<int>::iterator st_it = status.begin(); st_it != status.end(); ++st_it) {
0083           if (channelAllowed == *st_it || (channelAllowed < 0 && abs(channelAllowed) <= *st_it)) {
0084             passChannelLimitation = true;
0085             break;
0086           }
0087         }
0088       }
0089
0090       if (passChannelLimitation || limitDeadCellToChannelStatusEB.empty()) {
0091         if (bInfo.boundaryEnergy > highestEnergyDepositAroundDeadCell) {
0092           highestEnergyDepositAroundDeadCell = bInfo.boundaryET;
0093           //highestEnergyDepositAroundDeadCell = bInfo.boundaryEnergy;
0094         }
0095       }
0096     }
0097
0098     for (int i = 0; i < (int)v_boundaryInfoDeadCells_EE.size(); ++i) {
0099       BoundaryInformation bInfo = v_boundaryInfoDeadCells_EE[i];
0100
0101       //check if channel limitation rejectsbInfo
0102       bool passChannelLimitation = false;
0103       std::vector<int> status = bInfo.channelStatus;
0104
0105       for (int cs = 0; cs < (int)limitDeadCellToChannelStatusEE.size(); ++cs) {
0106         int channelAllowed = limitDeadCellToChannelStatusEE[cs];
0107
0108         for (std::vector<int>::iterator st_it = status.begin(); st_it != status.end(); ++st_it) {
0109           if (channelAllowed == *st_it || (channelAllowed < 0 && abs(channelAllowed) <= *st_it)) {
0110             passChannelLimitation = true;
0111             break;
0112           }
0113         }
0114       }
0115
0116       if (passChannelLimitation || limitDeadCellToChannelStatusEE.empty()) {
0117         if (bInfo.boundaryEnergy > highestEnergyDepositAroundDeadCell) {
0118           highestEnergyDepositAroundDeadCell = bInfo.boundaryET;
0119           //highestEnergyDepositAroundDeadCell = bInfo.boundaryEnergy;
0120         }
0121       }
0122     }
0123
0124     if (highestEnergyDepositAroundDeadCell > maxBoundaryEnergy) {
0125       //            cout << "<<<<<<<<<< List of EB  Boundary objects <<<<<<<<<<" << endl;
0126       //            for (int i = 0; i < (int) v_boundaryInfoDeadCells_EB.size(); ++i) {
0127       //               BoundaryInformation bInfo = v_boundaryInfoDeadCells_EB[i];
0128       //               cout << "no of neighbouring RecHits:" << bInfo.recHits.size() << endl;
0129       //               cout << "no of neighbouring DetIds:" << bInfo.detIds.size() << endl;
0130       //               cout << "boundary energy:" << bInfo.boundaryEnergy << endl;
0131       //               cout << "Channel stati: ";
0132       //               for (std::vector<int>::iterator it = bInfo.channelStatus.begin(); it != bInfo.channelStatus.end(); ++it) {
0133       //                  cout << *it << " ";
0134       //               }
0135       //               cout << endl;
0136       //            }
0137       //            cout << "<<<<<<<<<< List of EE  Boundary objects <<<<<<<<<<" << endl;
0138       //            for (int i = 0; i < (int) v_boundaryInfoDeadCells_EE.size(); ++i) {
0139       //               BoundaryInformation bInfo = v_boundaryInfoDeadCells_EE[i];
0140       //               cout << "no of neighbouring RecHits:" << bInfo.recHits.size() << endl;
0141       //               cout << "no of neighbouring DetIds:" << bInfo.detIds.size() << endl;
0142       //               cout << "boundary energy:" << bInfo.boundaryEnergy << endl;
0143       //               cout << "Channel stati: ";
0144       //               for (std::vector<int>::iterator it = bInfo.channelStatus.begin(); it != bInfo.channelStatus.end(); ++it) {
0145       //                  cout << *it << " ";
0146       //               }
0147       //               cout << endl;
0148       //            }
0149       return true;
0150     } else
0151       return false;
0152   }
0153
0154   bool isGapEcalCluster(double maxGapEnergyEB = 10, double maxGapEnergyEE = 10) const {
0155     float highestEnergyDepositAlongGapEB = 0;
0156
0157     for (int i = 0; i < (int)v_enNeighboursGap_EB.size(); ++i) {
0158       BoundaryInformation gapInfo = v_enNeighboursGap_EB[i];
0159
0160       if (gapInfo.boundaryEnergy > highestEnergyDepositAlongGapEB) {
0161         highestEnergyDepositAlongGapEB = gapInfo.boundaryET;
0162         //highestEnergyDepositAlongGapEB = gapInfo.boundaryEnergy;
0163       }
0164     }
0165
0166     float highestEnergyDepositAlongGapEE = 0;
0167
0168     for (int i = 0; i < (int)v_enNeighboursGap_EE.size(); ++i) {
0169       BoundaryInformation gapInfo = v_enNeighboursGap_EE[i];
0170
0171       if (gapInfo.boundaryEnergy > highestEnergyDepositAlongGapEE) {
0172         highestEnergyDepositAlongGapEE = gapInfo.boundaryET;
0173         //highestEnergyDepositAlongGapEE = gapInfo.boundaryEnergy;
0174       }
0175     }
0176
0177     if (highestEnergyDepositAlongGapEB > maxGapEnergyEB || highestEnergyDepositAlongGapEE > maxGapEnergyEE) {
0178       //            cout << "<<<<<<<<<< List of EB Gap objects <<<<<<<<<<" << endl;
0179       //            for (int i = 0; i < (int) v_enNeighboursGap_EB.size(); ++i) {
0180       //               BoundaryInformation gapInfo = v_enNeighboursGap_EB[i];
0181       //               cout << "no of neighbouring RecHits:" << gapInfo.recHits.size() << endl;
0182       //               cout << "no of neighbouring DetIds:" << gapInfo.detIds.size() << endl;
0183       //               cout << "gap energy:" << gapInfo.boundaryEnergy << endl;
0184       //            }
0185       //            cout << "<<<<<<<<<< List of EE Gap objects <<<<<<<<<<" << endl;
0186       //            for (int i = 0; i < (int) v_enNeighboursGap_EE.size(); ++i) {
0187       //               BoundaryInformation gapInfo = v_enNeighboursGap_EE[i];
0188       //               cout << "no of neighbouring RecHits:" << gapInfo.recHits.size() << endl;
0189       //               cout << "no of neighbouring DetIds:" << gapInfo.detIds.size() << endl;
0190       //               cout << "gap energy:" << gapInfo.boundaryEnergy << endl;
0191       //            }
0192       return true;
0193     } else
0194       return false;
0195   }
0196
0197   std::vector<BoundaryInformation> v_enNeighboursGap_EB;
0198   std::vector<BoundaryInformation> v_enNeighboursGap_EE;
0199
0200   std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EB;
0201   std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EE;
0202
0203 private:
0204 };
0205
0206 #endif /*ANOMALOUSECALVARIABLES_H_*/