AnomalousECALVariables.h

Go to the documentation of this file.

#ifndef ANOMALOUSECALVARIABLES_H_
#define ANOMALOUSECALVARIABLES_H_
//DataFormats/AnomalousEcalDataFormats/interface/AnomalousECALVariables.h
// system include files
#include <memory>

#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/METReco/interface/BoundaryInformation.h"

//using namespace edm;
//using namespace std;
/*
 * This class summarizes the information about the boundary energy calculated in EcalAnomalousEventFilter:
 * 1. next to ECAL border/gap
 * 2. next to masked ECAL channels: for each dead area with boundary energy above a threshold defined in the filter
 * the vector 'v_enDeadNeighbours_EB' or 'v_enDeadNeighbours_EE' is filled with the calculated boundary energy.
 * The determined size of the corresponding cluster is filled in v_enDeadNeighboursNoCells_EB/EE accordingly.
 *
 */
class AnomalousECALVariables {

   public:
      AnomalousECALVariables() {
         //energy next to ECAL Gap
         v_enNeighboursGap_EB.reserve(50);
         v_enNeighboursGap_EE.reserve(50);
         v_enNeighboursGap_EB.clear();
         v_enNeighboursGap_EE.clear();

         //energy around dead cells
         v_boundaryInfoDeadCells_EB = std::vector<BoundaryInformation> ();
         v_boundaryInfoDeadCells_EE = std::vector<BoundaryInformation> ();
         v_boundaryInfoDeadCells_EB.reserve(50);
         v_boundaryInfoDeadCells_EE.reserve(50);
         v_boundaryInfoDeadCells_EB.clear();
         v_boundaryInfoDeadCells_EE.clear();

      }
      ;

      AnomalousECALVariables(const std::vector<BoundaryInformation>& p_enNeighboursGap_EB,
            const std::vector<BoundaryInformation>& p_enNeighboursGap_EE, const std::vector<BoundaryInformation>& p_boundaryInfoDeadCells_EB,
            const std::vector<BoundaryInformation>& p_boundaryInfoDeadCells_EE) {

         v_boundaryInfoDeadCells_EB = std::vector<BoundaryInformation> ();
         v_boundaryInfoDeadCells_EE = std::vector<BoundaryInformation> ();
         v_boundaryInfoDeadCells_EB.reserve(50);
         v_boundaryInfoDeadCells_EE.reserve(50);
         v_boundaryInfoDeadCells_EB.clear();
         v_boundaryInfoDeadCells_EE.clear();
         v_boundaryInfoDeadCells_EB = p_boundaryInfoDeadCells_EB;
         v_boundaryInfoDeadCells_EE = p_boundaryInfoDeadCells_EE;

         v_enNeighboursGap_EB = p_enNeighboursGap_EB;
         v_enNeighboursGap_EE = p_enNeighboursGap_EE;
      }
      ;

      ~AnomalousECALVariables() {
         //cout<<"destructor AnomalousECAL"<<endl;
         v_enNeighboursGap_EB.clear();
         v_enNeighboursGap_EE.clear();
         v_boundaryInfoDeadCells_EB.clear();
         v_boundaryInfoDeadCells_EE.clear();
      }
      ;

      //returns true if at least 1 dead cell area was found in EcalAnomalousEventFilter with
      //boundary energy above threshold
      //Note: no sense to change this cut BELOW the threshold given in EcalAnomalousEventFilter

      bool isDeadEcalCluster(double maxBoundaryEnergy = 10,
            const std::vector<int>& limitDeadCellToChannelStatusEB = std::vector<int> (), const std::vector<int>& limitDeadCellToChannelStatusEE =
                  std::vector<int> ()) const {

         float highestEnergyDepositAroundDeadCell = 0;

         for (int i = 0; i < (int) v_boundaryInfoDeadCells_EB.size(); ++i) {
            BoundaryInformation bInfo = v_boundaryInfoDeadCells_EB[i];

            //check if channel limitation rejectsbInfo
            bool passChannelLimitation = false;
            std::vector<int> status = bInfo.channelStatus;

            for (int cs = 0; cs < (int) limitDeadCellToChannelStatusEB.size(); ++cs) {
               int channelAllowed = limitDeadCellToChannelStatusEB[cs];

               for (std::vector<int>::iterator st_it = status.begin(); st_it != status.end(); ++st_it) {

                  if (channelAllowed == *st_it || (channelAllowed < 0 && abs(channelAllowed) <= *st_it)) {
                     passChannelLimitation = true;
                     break;
                  }
               }
            }

            if (passChannelLimitation || limitDeadCellToChannelStatusEB.empty()) {

               if (bInfo.boundaryEnergy > highestEnergyDepositAroundDeadCell) {
                  highestEnergyDepositAroundDeadCell = bInfo.boundaryET;
                  //highestEnergyDepositAroundDeadCell = bInfo.boundaryEnergy;
               }
            }
         }

         for (int i = 0; i < (int) v_boundaryInfoDeadCells_EE.size(); ++i) {
            BoundaryInformation bInfo = v_boundaryInfoDeadCells_EE[i];

            //check if channel limitation rejectsbInfo
            bool passChannelLimitation = false;
            std::vector<int> status = bInfo.channelStatus;

            for (int cs = 0; cs < (int) limitDeadCellToChannelStatusEE.size(); ++cs) {
               int channelAllowed = limitDeadCellToChannelStatusEE[cs];

               for (std::vector<int>::iterator st_it = status.begin(); st_it != status.end(); ++st_it) {

                  if (channelAllowed == *st_it || (channelAllowed < 0 && abs(channelAllowed) <= *st_it)) {
                     passChannelLimitation = true;
                     break;
                  }
               }
            }

            if (passChannelLimitation || limitDeadCellToChannelStatusEE.empty()) {

               if (bInfo.boundaryEnergy > highestEnergyDepositAroundDeadCell){
                  highestEnergyDepositAroundDeadCell = bInfo.boundaryET;
                  //highestEnergyDepositAroundDeadCell = bInfo.boundaryEnergy;
               }
            }
         }

         if (highestEnergyDepositAroundDeadCell > maxBoundaryEnergy) {
            //            cout << "<<<<<<<<<< List of EB  Boundary objects <<<<<<<<<<" << endl;
            //            for (int i = 0; i < (int) v_boundaryInfoDeadCells_EB.size(); ++i) {
            //               BoundaryInformation bInfo = v_boundaryInfoDeadCells_EB[i];
            //               cout << "no of neighbouring RecHits:" << bInfo.recHits.size() << endl;
            //               cout << "no of neighbouring DetIds:" << bInfo.detIds.size() << endl;
            //               cout << "boundary energy:" << bInfo.boundaryEnergy << endl;
            //               cout << "Channel stati: ";
            //               for (std::vector<int>::iterator it = bInfo.channelStatus.begin(); it != bInfo.channelStatus.end(); ++it) {
            //                  cout << *it << " ";
            //               }
            //               cout << endl;
            //            }
            //            cout << "<<<<<<<<<< List of EE  Boundary objects <<<<<<<<<<" << endl;
            //            for (int i = 0; i < (int) v_boundaryInfoDeadCells_EE.size(); ++i) {
            //               BoundaryInformation bInfo = v_boundaryInfoDeadCells_EE[i];
            //               cout << "no of neighbouring RecHits:" << bInfo.recHits.size() << endl;
            //               cout << "no of neighbouring DetIds:" << bInfo.detIds.size() << endl;
            //               cout << "boundary energy:" << bInfo.boundaryEnergy << endl;
            //               cout << "Channel stati: ";
            //               for (std::vector<int>::iterator it = bInfo.channelStatus.begin(); it != bInfo.channelStatus.end(); ++it) {
            //                  cout << *it << " ";
            //               }
            //               cout << endl;
            //            }
            return true;
         } else
            return false;
      }

      bool isGapEcalCluster(double maxGapEnergyEB = 10, double maxGapEnergyEE = 10) const {

         float highestEnergyDepositAlongGapEB = 0;

         for (int i = 0; i < (int) v_enNeighboursGap_EB.size(); ++i) {
            BoundaryInformation gapInfo = v_enNeighboursGap_EB[i];

            if (gapInfo.boundaryEnergy > highestEnergyDepositAlongGapEB){
               highestEnergyDepositAlongGapEB = gapInfo.boundaryET;
               //highestEnergyDepositAlongGapEB = gapInfo.boundaryEnergy;
            }
         }

         float highestEnergyDepositAlongGapEE = 0;

         for (int i = 0; i < (int) v_enNeighboursGap_EE.size(); ++i) {
            BoundaryInformation gapInfo = v_enNeighboursGap_EE[i];

            if (gapInfo.boundaryEnergy > highestEnergyDepositAlongGapEE){
               highestEnergyDepositAlongGapEE = gapInfo.boundaryET;
               //highestEnergyDepositAlongGapEE = gapInfo.boundaryEnergy;
            }
         }

         if (highestEnergyDepositAlongGapEB > maxGapEnergyEB || highestEnergyDepositAlongGapEE > maxGapEnergyEE) {
            //            cout << "<<<<<<<<<< List of EB Gap objects <<<<<<<<<<" << endl;
            //            for (int i = 0; i < (int) v_enNeighboursGap_EB.size(); ++i) {
            //               BoundaryInformation gapInfo = v_enNeighboursGap_EB[i];
            //               cout << "no of neighbouring RecHits:" << gapInfo.recHits.size() << endl;
            //               cout << "no of neighbouring DetIds:" << gapInfo.detIds.size() << endl;
            //               cout << "gap energy:" << gapInfo.boundaryEnergy << endl;
            //            }
            //            cout << "<<<<<<<<<< List of EE Gap objects <<<<<<<<<<" << endl;
            //            for (int i = 0; i < (int) v_enNeighboursGap_EE.size(); ++i) {
            //               BoundaryInformation gapInfo = v_enNeighboursGap_EE[i];
            //               cout << "no of neighbouring RecHits:" << gapInfo.recHits.size() << endl;
            //               cout << "no of neighbouring DetIds:" << gapInfo.detIds.size() << endl;
            //               cout << "gap energy:" << gapInfo.boundaryEnergy << endl;
            //            }
            return true;
         } else
            return false;
      }

      std::vector<BoundaryInformation> v_enNeighboursGap_EB;
      std::vector<BoundaryInformation> v_enNeighboursGap_EE;

      std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EB;
      std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EE;

   private:

};

#endif /*ANOMALOUSECALVARIABLES_H_*/