#include <EcalGeomPhiSymHelper.h>

Public Member Functions
void	setup (const CaloGeometry geometry, const EcalChannelStatus chstatus, int statusThreshold)
Public Attributes
double	cellArea_ [kEndcWedgesX][kEndcWedgesY]
double	cellPhi_ [kEndcWedgesX][kEndcWedgesY]
GlobalPoint	cellPos_ [kEndcWedgesX][kEndcWedgesY]
int	endcapRing_ [kEndcWedgesX][kEndcWedgesY]
double	etaBoundary_ [kEndcEtaRings+1]
bool	goodCell_barl [kBarlRings][kBarlWedges][kSides]
bool	goodCell_endc [kEndcWedgesX][kEndcWedgesX][kSides]
double	meanCellArea_ [kEndcEtaRings]
int	nBads_barl [kBarlRings]
int	nBads_endc [kEndcEtaRings]
int	nRing_ [kEndcEtaRings]
double	phi_endc_ [kMaxEndciPhi][kEndcEtaRings]

Detailed Description

Definition at line 20 of file EcalGeomPhiSymHelper.h.

Member Function Documentation

void EcalGeomPhiSymHelper::setup	(	const CaloGeometry *	geometry,
		const EcalChannelStatus *	chstatus,
		int	statusThreshold
	)

Definition at line 19 of file EcalGeomPhiSymHelper.cc.

References abs, cellArea_, cellPhi_, cellPos_, gather_cfg::cout, DetId::Ecal, EcalBarrel, EcalEndcap, endcapRing_, eta(), PV3DBase< T, PVType, FrameType >::eta(), etaBoundary_, CaloCellGeometry::getCorners(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), CaloGeometry::getSubdetectorGeometry(), CaloGeometry::getValidDetIds(), goodCell_barl, goodCell_endc, EEDetId::hashedIndex(), i, EBDetId::ieta(), EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), kBarlRings, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kMaxEndciPhi, kSides, meanCellArea_, nBads_barl, nBads_endc, nRing_, dbtoconf::out, PV3DBase< T, PVType, FrameType >::phi(), phi, phi_endc_, relativeConstraints::ring, EEDetId::zside(), and EBDetId::zside().

Referenced by PhiSymmetryCalibration::setUp().

                                                    {

  
  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++) {

    const CaloCellGeometry *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;
              
    }
}