EcalCalibAlgos/src/EcalGeomPhiSymHelper.cc

0001 #include "Calibration/EcalCalibAlgos/interface/EcalGeomPhiSymHelper.h"
0002
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004
0005 // Geometry
0006 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0010 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0011
0012 //Channel status
0013
0014 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0015 #include "CondFormats/EcalObjects/interface/EcalChannelStatusCode.h"
0016 #include <fstream>
0017
0018 void EcalGeomPhiSymHelper::setup(const CaloGeometry* geometry, const EcalChannelStatus* chStatus, int statusThresold) {
0019   for (int ieta = 0; ieta < kBarlRings; ieta++)
0020     nBads_barl[ieta] = 0;
0021   for (int ring = 0; ring < kEndcEtaRings; ring++)
0022     nBads_endc[ring] = 0;
0023
0024   for (int ix = 0; ix < kEndcWedgesX; ix++) {
0025     for (int iy = 0; iy < kEndcWedgesY; iy++) {
0026       cellPhi_[ix][iy] = 0.;
0027       cellArea_[ix][iy] = 0.;
0028       endcapRing_[ix][iy] = -1;
0029     }
0030   }
0031
0032   // loop over all barrel crystals
0033   const std::vector<DetId>& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
0034   std::vector<DetId>::const_iterator barrelIt;
0035
0036   for (barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); barrelIt++) {
0037     EBDetId eb(*barrelIt);
0038
0039     int sign = eb.zside() > 0 ? 1 : 0;
0040
0041     int chs = (*chStatus)[*barrelIt].getStatusCode() & 0x001F;
0042     if (chs <= statusThresold)
0043       goodCell_barl[abs(eb.ieta()) - 1][eb.iphi() - 1][sign] = true;
0044
0045     if (!goodCell_barl[abs(eb.ieta()) - 1][eb.iphi() - 1][sign])
0046       nBads_barl[abs(eb.ieta()) - 1]++;
0047   }
0048
0049   const CaloSubdetectorGeometry* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0050
0051   for (int ix = 0; ix < kEndcWedgesX; ix++) {
0052     for (int iy = 0; iy < kEndcWedgesY; iy++) {
0053       cellPos_[ix][iy] = GlobalPoint(0., 0., 0.);
0054       cellPhi_[ix][iy] = 0.;
0055       cellArea_[ix][iy] = 0.;
0056       endcapRing_[ix][iy] = -1;
0057     }
0058   }
0059
0060   const std::vector<DetId>& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
0061   std::vector<DetId>::const_iterator endcapIt;
0062   for (endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); endcapIt++) {
0063     auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
0064     EEDetId ee(*endcapIt);
0065     int ix = ee.ix() - 1;
0066     int iy = ee.iy() - 1;
0067
0068     int sign = ee.zside() > 0 ? 1 : 0;
0069
0070     // store all crystal positions
0071     cellPos_[ix][iy] = cellGeometry->getPosition();
0072     cellPhi_[ix][iy] = cellGeometry->getPosition().phi();
0073
0074     // calculate and store eta-phi area for each crystal front face Shoelace formuls
0075     const CaloCellGeometry::CornersVec& cellCorners(cellGeometry->getCorners());
0076     cellArea_[ix][iy] = 0.;
0077
0078     for (int i = 0; i < 4; i++) {
0079       int iplus1 = i == 3 ? 0 : i + 1;
0080       cellArea_[ix][iy] += cellCorners[i].eta() * float(cellCorners[iplus1].phi()) -
0081                            cellCorners[iplus1].eta() * float(cellCorners[i].phi());
0082     }
0083
0084     cellArea_[ix][iy] = fabs(cellArea_[ix][iy]) / 2.;
0085     /*
0086     const double deltaPhi =
0087       (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaPhi(ee);
0088
0089     const double deltaEta =
0090       (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaEta(ee) ;
0091
0092     cellArea_[ix][iy] = deltaEta*deltaPhi;
0093 */
0094     int chs = (*chStatus)[*endcapIt].getStatusCode() & 0x001F;
0095     if (chs <= statusThresold)
0096       goodCell_endc[ix][iy][sign] = true;
0097   }
0098
0099   // get eta boundaries for each endcap ring
0100   etaBoundary_[0] = 1.479;
0101   etaBoundary_[39] = 3.;  //It was 4. !!!
0102   for (int ring = 1; ring < kEndcEtaRings; ring++) {
0103     double eta_ring_minus1 = cellPos_[ring - 1][50].eta();
0104     double eta_ring = cellPos_[ring][50].eta();
0105     etaBoundary_[ring] = (eta_ring + eta_ring_minus1) / 2.;
0106     std::cout << "Eta ring " << ring << " : " << eta_ring << std::endl;
0107   }
0108
0109   // determine to which ring each endcap crystal belongs,
0110   // the number of crystals in each ring,
0111   // and the mean eta-phi area of the crystals in each ring
0112
0113   for (int ring = 0; ring < kEndcEtaRings; ring++) {
0114     nRing_[ring] = 0;
0115     meanCellArea_[ring] = 0.;
0116     for (int ix = 0; ix < kEndcWedgesX; ix++) {
0117       for (int iy = 0; iy < kEndcWedgesY; iy++) {
0118         if (fabs(cellPos_[ix][iy].eta()) > etaBoundary_[ring] &&
0119             fabs(cellPos_[ix][iy].eta()) < etaBoundary_[ring + 1]) {
0120           meanCellArea_[ring] += cellArea_[ix][iy];
0121           endcapRing_[ix][iy] = ring;
0122           nRing_[ring]++;
0123
0124           for (int sign = 0; sign < kSides; sign++) {
0125             if (!goodCell_endc[ix][iy][sign])
0126               nBads_endc[ring]++;
0127           }  //sign
0128
0129         }  //if
0130       }  //ix
0131     }  //iy
0132
0133     meanCellArea_[ring] /= nRing_[ring];
0134
0135   }  //ring
0136
0137   // fill phi_endc[ip][ring] vector
0138   for (int ring = 0; ring < kEndcEtaRings; ring++) {
0139     for (int i = 0; i < kMaxEndciPhi; i++)
0140       phi_endc_[i][ring] = 0.;
0141
0142     float philast = -999.;
0143     for (int ip = 0; ip < nRing_[ring]; ip++) {
0144       float phimin = 999.;
0145       for (int ix = 0; ix < kEndcWedgesX; ix++) {
0146         for (int iy = 0; iy < kEndcWedgesY; iy++) {
0147           if (endcapRing_[ix][iy] == ring) {
0148             if (cellPhi_[ix][iy] < phimin && cellPhi_[ix][iy] > philast) {
0149               phimin = cellPhi_[ix][iy];
0150             }  //if edges
0151           }  //if ring
0152         }  //iy
0153       }  //ix
0154       phi_endc_[ip][ring] = phimin;
0155       philast = phimin;
0156     }  //ip
0157
0158   }  //ring
0159
0160   // Print out detid->ring association
0161   std::fstream eeringsf("endcaprings.dat", std::ios::out);
0162   for (endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); endcapIt++) {
0163     EEDetId eedet(*endcapIt);
0164     eeringsf << eedet.hashedIndex() << " " << endcapRing_[eedet.ix() - 1][eedet.iy() - 1] << " "
0165              << cellPhi_[eedet.ix() - 1][eedet.iy() - 1] << " "
0166              << cellArea_[eedet.ix() - 1][eedet.iy() - 1] / meanCellArea_[endcapRing_[eedet.ix() - 1][eedet.iy() - 1]]
0167              << std::endl;
0168   }
0169 }