PhotonFixCMS.cc

Go to the documentation of this file.

#include <cmath>
#include <cassert>
#include <fstream>
#include <iomanip>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/Records/interface/EcalBarrelGeometryRecord.h"
#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
#include "Geometry/Records/interface/EcalEndcapGeometryRecord.h"
#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"

#include "HiggsAnalysis/HiggsToGammaGamma/interface/PhotonFixCMS.h"

PhotonFixCMS::PhotonFixCMS(const reco::Photon &p):
  pf(p.energy(),p.superCluster()->eta(),p.superCluster()->phi(),p.r9()) {
}

bool PhotonFixCMS::initialise(const edm::EventSetup &iSetup, const std::string &s) {
  
  if(PhotonFix::initialised()) return false;
 
  PhotonFix::initialiseParameters(s);
  
  // Get ECAL geometry
  edm::ESHandle<CaloGeometry> geoHandle;
  iSetup.get<CaloGeometryRecord>().get(geoHandle);
  const CaloGeometry& geometry = *geoHandle;
  
  // EB
  const CaloSubdetectorGeometry *barrelGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
  
  double bc[170][360][2];
  for(int iz(0);iz<2;iz++) {
    for(int ie(0);ie<85;ie++) {
      int id = ie+1;
      if (iz==0) id = ie-85; 
      for(int ip(0);ip<360;ip++) {
    EBDetId eb(id,ip+1);
    
    auto cellGeometry = barrelGeometry->getGeometry(eb);
    const GlobalPoint& crystalPos = cellGeometry->getPosition();
    bc[85*iz+ie][ip][0]=crystalPos.eta();
    bc[85*iz+ie][ip][1]=crystalPos.phi();
      }
    }
  }
  
  for(unsigned i(0);i<169;i++) {
    for(unsigned j(0);j<360;j++) {
      unsigned k((j+1)%360);
     
      double eta = 0.25*(    bc[i][j][0]+bc[i+1][j][0]+
                     bc[i][k][0]+bc[i+1][k][0]);
      double phi = PhotonFix::GetaPhi(PhotonFix::GetaPhi(bc[i][j][1],bc[i+1][j][1]),
                PhotonFix::GetaPhi(bc[i][k][1],bc[i+1][k][1]));

      PhotonFix::barrelCGap(i,j,0,eta);
      PhotonFix::barrelCGap(i,j,1,phi);
      
      if((i%5)==4 && (j%2)==1) {
      PhotonFix::barrelSGap(i/5,j/2,0,eta);
    PhotonFix::barrelSGap(i/5,j/2,1,phi);   
      }
      
      if((j%20)==19) {
      if(i== 19) {PhotonFix::barrelMGap(0,j/20,0,eta); PhotonFix::barrelMGap(0,j/20,1,phi);}
      if(i== 39) {PhotonFix::barrelMGap(1,j/20,0,eta); PhotonFix::barrelMGap(1,j/20,1,phi);}
      if(i== 59) {PhotonFix::barrelMGap(2,j/20,0,eta); PhotonFix::barrelMGap(2,j/20,1,phi);}
      if(i== 84) {PhotonFix::barrelMGap(3,j/20,0,eta); PhotonFix::barrelMGap(3,j/20,1,phi);}
      if(i==109) {PhotonFix::barrelMGap(4,j/20,0,eta); PhotonFix::barrelMGap(4,j/20,1,phi);}
      if(i==129) {PhotonFix::barrelMGap(5,j/20,0,eta); PhotonFix::barrelMGap(5,j/20,1,phi);}
      if(i==149) {PhotonFix::barrelMGap(6,j/20,0,eta); PhotonFix::barrelMGap(6,j/20,1,phi);}
    
      }
    }
  }
  
  // EE
  const CaloSubdetectorGeometry *endcapGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
  
  double ec[2][100][100][2];
  bool valid[100][100];
  int val_count=0;
  for(int iz(0);iz<2;iz++) {
    for(int ix(0);ix<100;ix++) {
      for(int iy(0);iy<100;iy++) {
    valid[ix][iy] = EEDetId::validDetId(ix+1,iy+1,2*iz-1);
if(iz==0) PhotonFix::endcapCrystal(ix,iy,valid[ix][iy]);
    if(valid[ix][iy]) {
      EEDetId ee(ix+1,iy+1,2*iz-1);
    val_count+=1;
      
      auto cellGeometry = endcapGeometry->getGeometry(ee);
      const GlobalPoint& crystalPos = cellGeometry->getPosition();
      ec[iz][ix][iy][0]=asinh(crystalPos.x()/fabs(crystalPos.z()));
      ec[iz][ix][iy][1]=asinh(crystalPos.y()/fabs(crystalPos.z()));
    }
      }
    }
  }
  std::cout << "GG valid " << val_count << std::endl;
  double c[2];
  for(unsigned iz(0);iz<2;iz++) {
    unsigned nC(0),nS(0);
    for(unsigned i(0);i<99;i++) {
      for(unsigned j(0);j<99;j++) {
    if(valid[i][j  ] && valid[i+1][j  ] && 
       valid[i][j+1] && valid[i+1][j+1]) {
      for(unsigned k(0);k<2;k++) {

      c[k] = 0.25*(ec[iz][i][j][k]+ec[iz][i+1][j][k]+ec[iz][i][j+1][k]+ec[iz][i+1][j+1][k]);

        PhotonFix::endcapCGap(iz,nC,k,c[k]);     
    }
      
      if((i%5)==4 && (j%5)==4) {
        for(unsigned k(0);k<2;k++) {
        PhotonFix::endcapSGap(iz,nS,k,c[k]);     
        }
        nS++;
      }
      nC++;
    }
      }
    }
    std::cout << "Endcap number of crystal, submodule boundaries = "
          << nC << ", " << nS << std::endl;
  }
  
  // Hardcode EE D-module gap to 0,0
    PhotonFix::endcapMGap(0,0,0,0.0);    
    PhotonFix::endcapMGap(0,0,1,0.0);    
    PhotonFix::endcapMGap(1,0,0,0.0);    
    PhotonFix::endcapMGap(1,0,1,0.0);    
  
  return true;
}


double PhotonFixCMS::fixedEnergy() const {
  
  return pf.fixedEnergy();
}

double PhotonFixCMS::sigmaEnergy() const {
  
  return pf.sigmaEnergy();
}

const PhotonFix& PhotonFixCMS::photonFix() const {
  return pf;
}