HGCGeometryValidation.cc

Go to the documentation of this file.

#include <memory>
#include <iostream>

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/HcalSimNumberingRecord.h"
#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
#include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
#include "Geometry/HcalCommonData/interface/HcalCellType.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Utilities/interface/EDGetToken.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "DQMServices/Core/interface/DQMStore.h"

#include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
#include "SimDataFormats/ValidationFormats/interface/PHGCalValidInfo.h"
#include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"

#include "PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h"

#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Geometry/Vector3D.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"

const double mmtocm = 0.1;

class HGCGeometryValidation : public DQMEDAnalyzer {

public:

  explicit HGCGeometryValidation( const edm::ParameterSet& );
  ~HGCGeometryValidation() override;
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

protected:
  void dqmBeginRun(edm::Run const&, edm::EventSetup const&) override;
  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
  void analyze(const edm::Event&, const edm::EventSetup&) override;

private:
  edm::EDGetTokenT<PHGCalValidInfo> g4Token_;
  std::vector<std::string> geometrySource_;

  //HGCal geometry scheme
  std::vector<const HGCalDDDConstants*> hgcGeometry_;
  const HcalDDDSimConstants *hcons_;

  //histogram related stuff
  MonitorElement *heedzVsZ, *heedyVsY, *heedxVsX;
  MonitorElement *hefdzVsZ, *hefdyVsY, *hefdxVsX;
  MonitorElement *hebdzVsZ, *hebdyVsY, *hebdxVsX;
  MonitorElement *heedzVsLayer, *hefdzVsLayer, *hebdzVsLayer;
  MonitorElement *heedyVsLayer, *hefdyVsLayer, *hebdyVsLayer;
  MonitorElement *heedxVsLayer, *hefdxVsLayer, *hebdxVsLayer;
  MonitorElement *heeXG4VsId,   *hefXG4VsId,   *hebXG4VsId;
  MonitorElement *heeYG4VsId,   *hefYG4VsId,   *hebYG4VsId;
  MonitorElement *heeZG4VsId,   *hefZG4VsId,   *hebZG4VsId;
  MonitorElement *hebLayerVsEnStep, *hefLayerVsEnStep, *heeLayerVsEnStep;
    
  MonitorElement *heeTotEdepStep, *hefTotEdepStep, *hebTotEdepStep;
  MonitorElement *heedX, *heedY, *heedZ;
  MonitorElement *hefdX, *hefdY, *hefdZ;
  MonitorElement *hebdX, *hebdY, *hebdZ;
};

HGCGeometryValidation::HGCGeometryValidation(const edm::ParameterSet &cfg) : hcons_(nullptr) {

  g4Token_ = consumes<PHGCalValidInfo>(cfg.getParameter<edm::InputTag>("g4Source"));
  geometrySource_ = cfg.getUntrackedParameter< std::vector<std::string> >("geometrySource");
}

HGCGeometryValidation::~HGCGeometryValidation() { }

void HGCGeometryValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}

void HGCGeometryValidation::dqmBeginRun(const edm::Run& iRun, 
                    const edm::EventSetup& iSetup) {

  //initiating hgcnumbering
  for (size_t i=0; i<geometrySource_.size(); i++) {
    if (geometrySource_[i].find("Hcal") != std::string::npos) {
      edm::ESHandle<HcalDDDSimConstants> pHRNDC;
      iSetup.get<HcalSimNumberingRecord>().get(pHRNDC);
      if (pHRNDC.isValid()) {
    hcons_ = &(*pHRNDC);
    hgcGeometry_.push_back(nullptr);
      } else {
    edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for "
                      << geometrySource_[i];
      }
    } else {
      edm::ESHandle<HGCalDDDConstants> hgcGeom;
      iSetup.get<IdealGeometryRecord>().get(geometrySource_[i],hgcGeom);
      if (hgcGeom.isValid()) {
    hgcGeometry_.push_back(hgcGeom.product());
      } else {
    edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for "
                      << geometrySource_[i];
      }
    }
  }
}


void HGCGeometryValidation::bookHistograms(DQMStore::IBooker& iB, 
                                           edm::Run const&, 
                                           edm::EventSetup const&) {

  iB.setCurrentFolder("HGCAL/HGCalSimHitsV/Geometry");

  //initiating histograms
  heeTotEdepStep = iB.book1D("heeTotEdepStep","",100,0,100);
  hefTotEdepStep = iB.book1D("hefTotEdepStep","",100,0,100); 
  hebTotEdepStep = iB.book1D("hebTotEdepStep","",100,0,100);

  hebLayerVsEnStep = iB.book2D("hebLayerVsEnStep","",25,0,25,100,0,0.01);
  hefLayerVsEnStep = iB.book2D("hefLayerVsEnStep","",36,0,36,100,0,0.01);
  heeLayerVsEnStep = iB.book2D("heeLayerVsEnStep","",84,0,84,100,0,0.01);

  heeXG4VsId = iB.book2D("heeXG4VsId","",600,-300,300,600,-300,300);
  heeYG4VsId = iB.book2D("heeYG4VsId","",600,-300,300,600,-300,300);
  heeZG4VsId = iB.book2D("heeZG4VsId","",3000,320,350,3000,320,350);

  hefXG4VsId = iB.book2D("hefXG4VsId","",600,-300,300,600,-300,300);
  hefYG4VsId = iB.book2D("hefYG4VsId","",600,-300,300,600,-300,300);
  hefZG4VsId = iB.book2D("hefZG4VsId","",6000,350,410,6000,350,410);

  hebXG4VsId = iB.book2D("hebXG4VsId","",600,-300,300,600,-300,300);
  hebYG4VsId = iB.book2D("hebYG4VsId","",600,-300,300,600,-300,300);
  hebZG4VsId = iB.book2D("hebZG4VsId","",220,400,620,220,400,620);

  heedzVsZ = iB.book2D("heedzVsZ","",600,320,350,100,-1,1);
  heedyVsY = iB.book2D("heedyVsY","",400,-200,200,100,-1,1);
  heedxVsX = iB.book2D("heedxVsX","",400,-200,200,100,-1,1);

  hefdzVsZ = iB.book2D("hefdzVsZ","",1200,350,410,100,-1,1);
  hefdyVsY = iB.book2D("hefdyVsY","",400,-200,200,100,-1,1);
  hefdxVsX = iB.book2D("hefdxVsX","",400,-200,200,100,-1,1);

  hebdzVsZ = iB.book2D("hebdzVsZ","",220,400,620,100,-5,5);
  hebdyVsY = iB.book2D("hebdyVsY","",400,-200,200,100,-5,5);
  hebdxVsX = iB.book2D("hebdxVsX","",400,-200,200,100,-5,5);

  heedzVsLayer = iB.book2D("heedzVsLayer","",100,0,100,100,-1,1);
  hefdzVsLayer = iB.book2D("hefdzVsLayer","",40,0,40,100,-1,1);
  hebdzVsLayer = iB.book2D("hebdzVsLayer","",50,0,25,100,-5,5);

  heedyVsLayer = iB.book2D("heedyVsLayer","",100,0,100,100,-1,1);
  hefdyVsLayer = iB.book2D("hefdyVsLayer","",40,0,40,100,-1,1);
  hebdyVsLayer = iB.book2D("hebdyVsLayer","",50,0,25,100,-5,5);

  heedxVsLayer = iB.book2D("heedxVsLayer","",100,0,100,100,-1,1);
  hefdxVsLayer = iB.book2D("hefdxVsLayer","",40,0,40,500,-1,1);
  hebdxVsLayer = iB.book2D("hebdxVsLayer","",50,0,25,500,-5,5.0);

  heedX = iB.book1D("heedX","",100,-1,1); 
  heedY = iB.book1D("heedY","",100,-1,1);
  heedZ = iB.book1D("heedZ","",100,-1,1);

  hefdX = iB.book1D("hefdX","",100,-1,1); 
  hefdY = iB.book1D("hefdY","",100,-1,1);
  hefdZ = iB.book1D("hefdZ","",100,-1,1);

  hebdX = iB.book1D("hebdX","",100,-1,1); 
  hebdY = iB.book1D("hebdY","",100,-1,1);
  hebdZ = iB.book1D("hebdZ","",100,-1,1);
}

void HGCGeometryValidation::analyze(const edm::Event &iEvent, 
                    const edm::EventSetup &iSetup) {

  //Accessing G4 information
  edm::Handle<PHGCalValidInfo> infoLayer;
  iEvent.getByToken(g4Token_,infoLayer);

  if (infoLayer.isValid()) {
    //step vertex information
    std::vector<float> hitVtxX = infoLayer->hitvtxX();
    std::vector<float> hitVtxY = infoLayer->hitvtxY();
    std::vector<float> hitVtxZ = infoLayer->hitvtxZ();
    std::vector<unsigned int> hitDet = infoLayer->hitDets();
    std::vector<unsigned int> hitIdx = infoLayer->hitIndex();

    //energy information
    std::vector<float> edepLayerEE = infoLayer->eehgcEdep();
    std::vector<float> edepLayerHE = infoLayer->hefhgcEdep();
    std::vector<float> edepLayerHB = infoLayer->hebhgcEdep();
    
    unsigned int i;     
    for(i=0; i<edepLayerEE.size(); i++) {
      heeLayerVsEnStep->Fill(i,edepLayerEE.at(i));
    }
    
    for(i=0; i<edepLayerHE.size(); i++) {
      hefLayerVsEnStep->Fill(i,edepLayerHE.at(i));
    }

    for(i=0; i<edepLayerHB.size(); i++) {
      hebLayerVsEnStep->Fill(i,edepLayerHB.at(i));
    }

    //fill total energy deposited
    heeTotEdepStep->Fill((double)infoLayer->eeTotEdep());
    hefTotEdepStep->Fill((double)infoLayer->hefTotEdep()); 
    hebTotEdepStep->Fill((double)infoLayer->hebTotEdep());
    
    //loop over all hits
    for (unsigned int i=0; i<hitVtxX.size(); i++) {

      hitVtxX.at(i) = mmtocm*hitVtxX.at(i);
      hitVtxY.at(i) = mmtocm*hitVtxY.at(i);
      hitVtxZ.at(i) = mmtocm*hitVtxZ.at(i);

      if ((hitDet.at(i) == (unsigned int)(DetId::Forward)) ||
      (hitDet.at(i) == (unsigned int)(DetId::HGCalEE)) ||
      (hitDet.at(i) == (unsigned int)(DetId::HGCalHSi)) ||
      (hitDet.at(i) == (unsigned int)(DetId::HGCalHSc))) {
    double xx, yy;
    int    dtype(0), layer(0), zside(1);
    std::pair<float,float> xy;
    if (hitDet.at(i) == (unsigned int)(DetId::Forward)) {
      int subdet, wafer, celltype, cell;
      HGCalTestNumbering::unpackHexagonIndex(hitIdx.at(i), subdet, zside, 
                         layer, wafer, celltype,cell);
      dtype = (subdet==(int)(HGCEE)) ? 0 : 1; 
      xy = hgcGeometry_[dtype]->locateCell(cell,layer,wafer,true); //cm
    } else if ((hitDet.at(i) == (unsigned int)(DetId::HGCalEE)) ||
           (hitDet.at(i) == (unsigned int)(DetId::HGCalHSi))) {
      HGCSiliconDetId id(hitIdx.at(i));
      dtype = (id.det() == DetId::HGCalEE) ? 0 : 1;
      layer = id.layer();
      zside = id.zside();
      xy    = hgcGeometry_[dtype]->locateCell(layer,id.waferU(),
                          id.waferV(),id.cellU(),
                          id.cellV(),true,true);
    } else {
      HGCScintillatorDetId id(hitIdx.at(i));
      dtype = 2;
      layer = id.layer();
      zside = id.zside();
      xy    = hgcGeometry_[dtype]->locateCellTrap(layer,id.ietaAbs(),
                              id.iphi(), true);
    }
    double zz = hgcGeometry_[dtype]->waferZ(layer,true); //cm 
    if (zside < 0) zz = -zz;
    xx = (zside<0) ? -xy.first : xy.first; 
    yy = xy.second;

    if (dtype == 0) {

      heedzVsZ->Fill(zz, (hitVtxZ.at(i)-zz));
      heedyVsY->Fill(yy, (hitVtxY.at(i)-yy));
      heedxVsX->Fill(xx, (hitVtxX.at(i)-xx));
      
      heeXG4VsId->Fill(hitVtxX.at(i),xx);
      heeYG4VsId->Fill(hitVtxY.at(i),yy);
      heeZG4VsId->Fill(hitVtxZ.at(i),zz);

      heedzVsLayer->Fill(layer,(hitVtxZ.at(i)-zz));
      heedyVsLayer->Fill(layer,(hitVtxY.at(i)-yy));
      heedxVsLayer->Fill(layer,(hitVtxX.at(i)-xx));
      
      heedX->Fill((hitVtxX.at(i)-xx));
      heedZ->Fill((hitVtxZ.at(i)-zz));
      heedY->Fill((hitVtxY.at(i)-yy));

    } else if (dtype == 1) {
      
      hefdzVsZ->Fill(zz, (hitVtxZ.at(i)-zz));
      hefdyVsY->Fill(yy, (hitVtxY.at(i)-yy));
      hefdxVsX->Fill(xx, (hitVtxX.at(i)-xx));
      
      hefXG4VsId->Fill(hitVtxX.at(i),xx);
      hefYG4VsId->Fill(hitVtxY.at(i),yy);
      hefZG4VsId->Fill(hitVtxZ.at(i),zz);
      
      hefdzVsLayer->Fill(layer,(hitVtxZ.at(i)-zz));
      hefdyVsLayer->Fill(layer,(hitVtxY.at(i)-yy));
      hefdxVsLayer->Fill(layer,(hitVtxX.at(i)-xx));
      
      hefdX->Fill((hitVtxX.at(i)-xx));
      hefdZ->Fill((hitVtxZ.at(i)-zz));
      hefdY->Fill((hitVtxY.at(i)-yy));
    
    } else {

      hebdzVsZ->Fill(zz, (hitVtxZ.at(i)-zz));
      hebdyVsY->Fill(yy, (hitVtxY.at(i)-yy));
      hebdxVsX->Fill(xx, (hitVtxX.at(i)-xx));

      hebXG4VsId->Fill(hitVtxX.at(i),xx);
      hebYG4VsId->Fill(hitVtxY.at(i),yy);
      hebZG4VsId->Fill(hitVtxZ.at(i),zz);

      hebdzVsLayer->Fill(layer,(hitVtxZ.at(i)-zz));
      hebdyVsLayer->Fill(layer,(hitVtxY.at(i)-yy));
      hebdxVsLayer->Fill(layer,(hitVtxX.at(i)-xx));

      hebdX->Fill((hitVtxX.at(i)-xx));
      hebdZ->Fill((hitVtxZ.at(i)-zz));
      hebdY->Fill((hitVtxY.at(i)-yy));

    }
    
      } else if (hitDet.at(i) == (unsigned int)(DetId::Hcal)) {

    int subdet, zside, depth, eta, phi, lay;
    HcalTestNumbering::unpackHcalIndex(hitIdx.at(i), subdet, zside, depth, eta, phi, lay);
    HcalCellType::HcalCell cell = hcons_->cell(subdet, zside, lay, eta, phi);
    
    double zz = mmtocm*cell.rz; //mm --> cm
    if (zside == 0) zz = -zz;
    double rho = zz*tan(2.0*atan(exp(-cell.eta)));
    double xx  = rho * cos(cell.phi); //cm
    double yy  = rho * sin(cell.phi); //cm

    hebdzVsZ->Fill(zz, (hitVtxZ.at(i)-zz));
    hebdyVsY->Fill(yy, (hitVtxY.at(i)-yy));
    hebdxVsX->Fill(xx, (hitVtxX.at(i)-xx));

    hebXG4VsId->Fill(hitVtxX.at(i),xx);
    hebYG4VsId->Fill(hitVtxY.at(i),yy);
    hebZG4VsId->Fill(hitVtxZ.at(i),zz);

    hebdzVsLayer->Fill(lay,(hitVtxZ.at(i)-zz));
    hebdyVsLayer->Fill(lay,(hitVtxY.at(i)-yy));
    hebdxVsLayer->Fill(lay,(hitVtxX.at(i)-xx));

    hebdX->Fill((hitVtxX.at(i)-xx));
    hebdZ->Fill((hitVtxZ.at(i)-zz));
    hebdY->Fill((hitVtxY.at(i)-yy));
      }

    }//end G4 hits
    
  } else {
    edm::LogVerbatim("HGCalValid") << "HGCGeometryValidation::No PHGCalInfo ";
  }

}

//define this as a plug-in
DEFINE_FWK_MODULE(HGCGeometryValidation);