EcalGeomPhiSymHelper.cc

Go to the documentation of this file.

#include "Calibration/EcalCalibAlgos/interface/EcalGeomPhiSymHelper.h"

#include "FWCore/Framework/interface/ESHandle.h"

// Geometry
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"

//Channel status

#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
#include "CondFormats/EcalObjects/interface/EcalChannelStatusCode.h"
#include <fstream>

void EcalGeomPhiSymHelper::setup(const CaloGeometry* geometry, const EcalChannelStatus* chStatus, int statusThresold) {
  for (int ieta = 0; ieta < kBarlRings; ieta++)
    nBads_barl[ieta] = 0;
  for (int ring = 0; ring < kEndcEtaRings; ring++)
    nBads_endc[ring] = 0;

  for (int ix = 0; ix < kEndcWedgesX; ix++) {
    for (int iy = 0; iy < kEndcWedgesY; iy++) {
      cellPhi_[ix][iy] = 0.;
      cellArea_[ix][iy] = 0.;
      endcapRing_[ix][iy] = -1;
    }
  }

  // loop over all barrel crystals
  const std::vector<DetId>& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
  std::vector<DetId>::const_iterator barrelIt;

  for (barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); barrelIt++) {
    EBDetId eb(*barrelIt);

    int sign = eb.zside() > 0 ? 1 : 0;

    int chs = (*chStatus)[*barrelIt].getStatusCode() & 0x001F;
    if (chs <= statusThresold)
      goodCell_barl[abs(eb.ieta()) - 1][eb.iphi() - 1][sign] = true;

    if (!goodCell_barl[abs(eb.ieta()) - 1][eb.iphi() - 1][sign])
      nBads_barl[abs(eb.ieta()) - 1]++;
  }

  const CaloSubdetectorGeometry* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);

  for (int ix = 0; ix < kEndcWedgesX; ix++) {
    for (int iy = 0; iy < kEndcWedgesY; iy++) {
      cellPos_[ix][iy] = GlobalPoint(0., 0., 0.);
      cellPhi_[ix][iy] = 0.;
      cellArea_[ix][iy] = 0.;
      endcapRing_[ix][iy] = -1;
    }
  }

  const std::vector<DetId>& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
  std::vector<DetId>::const_iterator endcapIt;
  for (endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); endcapIt++) {
    auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
    EEDetId ee(*endcapIt);
    int ix = ee.ix() - 1;
    int iy = ee.iy() - 1;

    int sign = ee.zside() > 0 ? 1 : 0;

    // store all crystal positions
    cellPos_[ix][iy] = cellGeometry->getPosition();
    cellPhi_[ix][iy] = cellGeometry->getPosition().phi();

    // calculate and store eta-phi area for each crystal front face Shoelace formuls
    const CaloCellGeometry::CornersVec& cellCorners(cellGeometry->getCorners());
    cellArea_[ix][iy] = 0.;

    for (int i = 0; i < 4; i++) {
      int iplus1 = i == 3 ? 0 : i + 1;
      cellArea_[ix][iy] += cellCorners[i].eta() * float(cellCorners[iplus1].phi()) -
                           cellCorners[iplus1].eta() * float(cellCorners[i].phi());
    }

    cellArea_[ix][iy] = fabs(cellArea_[ix][iy]) / 2.;
    /*
    const double deltaPhi =
      (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaPhi(ee);

    const double deltaEta =
      (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaEta(ee) ;

    cellArea_[ix][iy] = deltaEta*deltaPhi;
*/
    int chs = (*chStatus)[*endcapIt].getStatusCode() & 0x001F;
    if (chs <= statusThresold)
      goodCell_endc[ix][iy][sign] = true;
  }

  // get eta boundaries for each endcap ring
  etaBoundary_[0] = 1.479;
  etaBoundary_[39] = 3.;  //It was 4. !!!
  for (int ring = 1; ring < kEndcEtaRings; ring++) {
    double eta_ring_minus1 = cellPos_[ring - 1][50].eta();
    double eta_ring = cellPos_[ring][50].eta();
    etaBoundary_[ring] = (eta_ring + eta_ring_minus1) / 2.;
    std::cout << "Eta ring " << ring << " : " << eta_ring << std::endl;
  }

  // determine to which ring each endcap crystal belongs,
  // the number of crystals in each ring,
  // and the mean eta-phi area of the crystals in each ring

  for (int ring = 0; ring < kEndcEtaRings; ring++) {
    nRing_[ring] = 0;
    meanCellArea_[ring] = 0.;
    for (int ix = 0; ix < kEndcWedgesX; ix++) {
      for (int iy = 0; iy < kEndcWedgesY; iy++) {
        if (fabs(cellPos_[ix][iy].eta()) > etaBoundary_[ring] &&
            fabs(cellPos_[ix][iy].eta()) < etaBoundary_[ring + 1]) {
          meanCellArea_[ring] += cellArea_[ix][iy];
          endcapRing_[ix][iy] = ring;
          nRing_[ring]++;

          for (int sign = 0; sign < kSides; sign++) {
            if (!goodCell_endc[ix][iy][sign])
              nBads_endc[ring]++;
          }  //sign

        }  //if
      }    //ix
    }      //iy

    meanCellArea_[ring] /= nRing_[ring];

  }  //ring

  // fill phi_endc[ip][ring] vector
  for (int ring = 0; ring < kEndcEtaRings; ring++) {
    for (int i = 0; i < kMaxEndciPhi; i++)
      phi_endc_[i][ring] = 0.;

    float philast = -999.;
    for (int ip = 0; ip < nRing_[ring]; ip++) {
      float phimin = 999.;
      for (int ix = 0; ix < kEndcWedgesX; ix++) {
        for (int iy = 0; iy < kEndcWedgesY; iy++) {
          if (endcapRing_[ix][iy] == ring) {
            if (cellPhi_[ix][iy] < phimin && cellPhi_[ix][iy] > philast) {
              phimin = cellPhi_[ix][iy];
            }  //if edges
          }    //if ring
        }      //iy
      }        //ix
      phi_endc_[ip][ring] = phimin;
      philast = phimin;
    }  //ip

  }  //ring

  // Print out detid->ring association
  std::fstream eeringsf("endcaprings.dat", std::ios::out);
  for (endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); endcapIt++) {
    EEDetId eedet(*endcapIt);
    eeringsf << eedet.hashedIndex() << " " << endcapRing_[eedet.ix() - 1][eedet.iy() - 1] << " "
             << cellPhi_[eedet.ix() - 1][eedet.iy() - 1] << " "
             << cellArea_[eedet.ix() - 1][eedet.iy() - 1] / meanCellArea_[endcapRing_[eedet.ix() - 1][eedet.iy() - 1]]
             << std::endl;
  }
}