HcalHBHEMuonSimAnalyzer.cc

Go to the documentation of this file.

#include <memory>
#include <iostream>
#include <vector>

#include <TTree.h>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"

#include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
#include "Calibration/IsolatedParticles/interface/eECALMatrix.h" 
#include "Calibration/IsolatedParticles/interface/eHCALMatrix.h" 

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#define EDM_ML_DEBUG

class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::SharedResources> {

public:
  explicit HcalHBHEMuonSimAnalyzer(const edm::ParameterSet&);
  ~HcalHBHEMuonSimAnalyzer() override;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  void beginJob() override;
  void analyze(edm::Event const&, edm::EventSetup const&) override;
  void beginRun(edm::Run const&, edm::EventSetup const&) override;
  void endRun(edm::Run const&, edm::EventSetup const&) override {}
  virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
  virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
  void         clearVectors();
  unsigned int matchId(const HcalDetId&, const HcalDetId&);
  double       activeLength(const DetId&);
  
  std::string                                   g4Label_, ebLabel_, eeLabel_;
  std::string                                   hcLabel_;
  int                                           verbosity_, maxDepth_;
  double                                        etaMax_;
  std::vector<HcalDDDRecConstants::HcalActiveLength> actHB_, actHE_;

  double                                        tMinE_, tMaxE_, tMinH_, tMaxH_;
  edm::Service<TFileService>                    fs_;
  edm::EDGetTokenT<edm::SimTrackContainer>      tok_SimTk_;
  edm::EDGetTokenT<edm::SimVertexContainer>     tok_SimVtx_;
  edm::EDGetTokenT<edm::PCaloHitContainer>      tok_caloEB_, tok_caloEE_;
  edm::EDGetTokenT<edm::PCaloHitContainer>      tok_caloHH_;

  static const int          depthMax_=7;
  const int                 idMuon_=13;
  TTree                    *tree_;
  unsigned int              runNumber_, eventNumber_, lumiNumber_, bxNumber_;
  std::vector<double>       ptGlob_, etaGlob_, phiGlob_, pMuon_;
  std::vector<double>       ecal3x3Energy_, hcal1x1Energy_;
  std::vector<unsigned int> ecalDetId_,     hcalDetId_,  hcalHot_;
  std::vector<double>       matchedId_;
  std::vector<double>       hcalDepthEnergy_[depthMax_];
  std::vector<double>       hcalDepthActiveLength_[depthMax_];
  std::vector<double>       hcalDepthEnergyHot_[depthMax_];
  std::vector<double>       hcalDepthActiveLengthHot_[depthMax_];
  std::vector<double>       hcalActiveLength_, hcalActiveLengthHot_;
};

HcalHBHEMuonSimAnalyzer::HcalHBHEMuonSimAnalyzer(const edm::ParameterSet& iConfig) {

  usesResource(TFileService::kSharedResource);

  //now do what ever initialization is needed
  g4Label_   = iConfig.getParameter<std::string>("ModuleLabel");
  ebLabel_   = iConfig.getParameter<std::string>("EBCollection");
  eeLabel_   = iConfig.getParameter<std::string>("EECollection");
  hcLabel_   = iConfig.getParameter<std::string>("HCCollection");
  verbosity_ = iConfig.getUntrackedParameter<int>("Verbosity",0);
  maxDepth_  = iConfig.getUntrackedParameter<int>("MaxDepth",4);
  etaMax_    = iConfig.getUntrackedParameter<double>("EtaMax", 3.0);
  tMinE_     = iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.);
  tMaxE_     = iConfig.getUntrackedParameter<double>("TimeMaxCutECAL",  500.);
  tMinH_     = iConfig.getUntrackedParameter<double>("TimeMinCutHCAL", -500.);
  tMaxH_     = iConfig.getUntrackedParameter<double>("TimeMaxCutHCAL",  500.);

  tok_SimTk_  = consumes<edm::SimTrackContainer>(edm::InputTag(g4Label_));
  tok_SimVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag(g4Label_));
  tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_,ebLabel_));
  tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_,eeLabel_));
  tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_,hcLabel_));
  if      (maxDepth_ > depthMax_) maxDepth_ = depthMax_;
  else if (maxDepth_ < 1)         maxDepth_ = 4;
#ifdef EDM_ML_DEBUG
  std::cout << "Labels: " << g4Label_ << ":" << ebLabel_ << ":" << eeLabel_
        << ":" << hcLabel_ << "\nVerbosity " << verbosity_ << " MaxDepth "
        << maxDepth_ << " Maximum Eta " << etaMax_ << " tMin|tMax "
        << tMinE_ << ":" << tMaxE_ << ":" << tMinH_ << ":" << tMaxH_
        << std::endl;
#endif
}

HcalHBHEMuonSimAnalyzer::~HcalHBHEMuonSimAnalyzer() { }

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

  clearVectors();
  bool debug(false);
#ifdef EDM_ML_DEBUG
  debug = ((verbosity_/10)>0);
#endif

  edm::ESHandle<HcalDDDRecConstants> pHRNDC;
  iSetup.get<HcalRecNumberingRecord>().get(pHRNDC);
  const HcalDDDRecConstants* hcons = &(*pHRNDC);

  runNumber_   = iEvent.id().run();
  eventNumber_ = iEvent.id().event();
  lumiNumber_  = iEvent.id().luminosityBlock();
  bxNumber_    = iEvent.bunchCrossing();
  
  //get Handles to SimTracks and SimHits
  edm::Handle<edm::SimTrackContainer> SimTk;
  iEvent.getByToken(tok_SimTk_,SimTk);
  edm::SimTrackContainer::const_iterator simTrkItr;
  edm::Handle<edm::SimVertexContainer> SimVtx;
  iEvent.getByToken(tok_SimVtx_,SimVtx);

  //get Handles to PCaloHitContainers of eb/ee/hbhe
  edm::Handle<edm::PCaloHitContainer> pcaloeb;
  iEvent.getByToken(tok_caloEB_, pcaloeb);
  edm::Handle<edm::PCaloHitContainer> pcaloee;
  iEvent.getByToken(tok_caloEE_, pcaloee);
  edm::Handle<edm::PCaloHitContainer> pcalohh;
  iEvent.getByToken(tok_caloHH_, pcalohh);
  std::vector<PCaloHit> calohh;
  bool testN(false);
  for (unsigned int k=1; k<pcalohh->size(); ++k) {
    // if it is a standard DetId bits 28..31 will carry the det #
    // for HCAL det # is 4 and if there is at least one hit in the collection
    // have det # which is not 4 this collection is created using TestNumbering
    int det = ((((*pcalohh)[k].id())>>28)&0xF);
    if (det != 4) {testN = true; break;}
  }
  if (testN) {
    for (edm::PCaloHitContainer::const_iterator itr=pcalohh->begin(); itr != pcalohh->end(); ++itr) {
      PCaloHit hit(*itr);
      DetId newid = HcalHitRelabeller::relabel(hit.id(),hcons);
#ifdef EDM_ML_DEBUG
      std::cout << "Old ID " << std::hex << hit.id() << std::dec << " New " 
        << HcalDetId(newid) << std::endl;
#endif
      hit.setID(newid.rawId());
      calohh.push_back(hit);
    }
  } else {
    calohh.insert(calohh.end(),pcalohh->begin(),pcalohh->end());
  }
  
  // get handles to calogeometry and calotopology
  edm::ESHandle<CaloGeometry> pG;
  iSetup.get<CaloGeometryRecord>().get(pG);
  const CaloGeometry* geo = pG.product();

  edm::ESHandle<MagneticField> bFieldH;
  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
  const MagneticField* bField = bFieldH.product();

  edm::ESHandle<CaloTopology> theCaloTopology;
  iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
  const CaloTopology *caloTopology = theCaloTopology.product();

  edm::ESHandle<HcalTopology> htopo;
  iSetup.get<HcalRecNumberingRecord>().get(htopo);
  const HcalTopology* theHBHETopology = htopo.product();

  // Loop over all SimTracks
  for (edm::SimTrackContainer::const_iterator simTrkItr=SimTk->begin();
       simTrkItr!= SimTk->end(); simTrkItr++) {
    if ((std::abs(simTrkItr->type()) == idMuon_) && (simTrkItr->vertIndex() == 0) &&
    (std::abs(simTrkItr->momentum().eta()) < etaMax_)) {
      unsigned int thisTrk = simTrkItr->trackId();
      spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, debug);

      double eEcal(0), eHcal(0), activeLengthTot(0), activeLengthHotTot(0);
      double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
      double activeL[depthMax_], activeHotL[depthMax_];
      unsigned int isHot(0);
      bool         tmpmatch(false); 
      for (int i=0; i<depthMax_; ++i) 
    eHcalDepth[i] = eHcalDepthHot[i] = activeL[i] = activeHotL[i] = -10000;

#ifdef EDM_ML_DEBUG
      if ((verbosity_%10) > 0)
    std::cout << "Track Type " << simTrkItr->type() << " Vertex "
          << simTrkItr->vertIndex() << " Charge " << simTrkItr->charge()
          << " Momentum " << simTrkItr->momentum().P() << ":"
          << simTrkItr->momentum().eta() << ":"
          << simTrkItr->momentum().phi() << " ECAL|HCAL " << trkD.okECAL
          << ":" << trkD.okHCAL << " Point " << trkD.pointECAL << ":"
          << trkD.pointHCAL << " Direction " << trkD.directionECAL.eta()
          << ":" << trkD.directionECAL.phi() << " | "
          << trkD.directionHCAL.eta() << ":" << trkD.directionHCAL.phi()
          << std::endl;
#endif
      bool propageback(false);
      spr::propagatedTrackDirection trkD_back = spr::propagateHCALBack(thisTrk, SimTk, SimVtx, geo, bField, debug);
      HcalDetId closestCell_back;
      if(trkD_back.okHCAL) {
    closestCell_back = (HcalDetId) (trkD_back.detIdHCAL);
    propageback = true;
      }
      if (trkD.okHCAL) {
    // Muon properties
    spr::trackAtOrigin tkvx = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
    ptGlob_.push_back(tkvx.momentum.perp());
    etaGlob_.push_back(tkvx.momentum.eta());
    phiGlob_.push_back(tkvx.momentum.phi());
    pMuon_.push_back(tkvx.momentum.mag());
#ifdef EDM_ML_DEBUG
    if ((verbosity_%10) > 0)
      std::cout << "Track at vertex " << tkvx.ok << " position " 
            << tkvx.position << " Momentum " << tkvx.momentum.mag()
            << ":" << tkvx.momentum.eta() << ":"
            << tkvx.momentum.phi() << " Charge " << tkvx.charge
            << std::endl;
#endif

    
    // Energy in ECAL
    DetId isoCell;
    if (trkD.okECAL) {
      isoCell = trkD.detIdECAL;
      eEcal = spr::eECALmatrix(isoCell, pcaloeb, pcaloee, geo, caloTopology, 1, 1, -100.0, -100.0, tMinE_, tMaxE_, debug);
    }

    // Energy in  Hcal
    const DetId closestCell(trkD.detIdHCAL);
    if ((propageback) && 
        (HcalDetId(closestCell).ieta()==HcalDetId(closestCell_back).ieta()) && 
        (HcalDetId(closestCell).iphi() == HcalDetId(closestCell_back).iphi()))
      tmpmatch= true;
    
    eHcal = spr::eHCALmatrix(theHBHETopology, closestCell, calohh,0,0, false, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, debug);
#ifdef EDM_ML_DEBUG
    if ((verbosity_%10) > 0)
      std::cout << "eEcal " << trkD.okECAL << ":" << eEcal << " eHcal " 
            << eHcal << std::endl;
#endif

    HcalSubdetector subdet = HcalDetId(closestCell).subdet();
    int             ieta   = HcalDetId(closestCell).ieta();
    int             iphi   = HcalDetId(closestCell).iphi();
    int             zside  = HcalDetId(closestCell).zside();
    bool            hbhe   = (std::abs(ieta) == 16);
    int             depthHE= hcons->getMinDepth(1,16,iphi,zside);
    std::vector<std::pair<double,int> > ehdepth;
    spr::energyHCALCell((HcalDetId)closestCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, -500.0, 500.0, depthHE, debug);
    for (unsigned int i=0; i<ehdepth.size(); ++i) {
      eHcalDepth[ehdepth[i].second-1] = ehdepth[i].first;
      HcalSubdetector subdet0 = (hbhe) ? ((ehdepth[i].second >= depthHE) ? HcalEndcap : HcalBarrel) : subdet;
      HcalDetId hcid0(subdet0,ieta,iphi,ehdepth[i].second);
      double actL = activeLength(DetId(hcid0));
      activeL[ehdepth[i].second-1] = actL;
      activeLengthTot += actL;
#ifdef EDM_ML_DEBUG
      if ((verbosity_%10) > 0)
        std::cout << hcid0 << " E " << ehdepth[i].first << " L " << actL 
              << std::endl;
#endif
    }

    HcalDetId hotCell;
#ifdef EDM_ML_DEBUG
    double h3x3 = 
#endif
      spr::eHCALmatrix(geo,theHBHETopology, closestCell, calohh, 1,1, hotCell, debug);
    isHot = matchId(closestCell,hotCell);
#ifdef EDM_ML_DEBUG
    if ((verbosity_%10) > 0)
      std::cout << "hcal 3X3  < " << h3x3 << ">" << " ClosestCell <" 
            << (HcalDetId)(closestCell) << "> hotCell id < " << hotCell
            << "> isHot" << isHot << std::endl;
#endif

    if (hotCell != HcalDetId()) {
      subdet = HcalDetId(hotCell).subdet();
      ieta   = HcalDetId(hotCell).ieta();
      iphi   = HcalDetId(hotCell).iphi();
      zside  = HcalDetId(hotCell).zside();
      hbhe   = (std::abs(ieta) == 16);
      depthHE= hcons->getMinDepth(1,16,iphi,zside);
      std::vector<std::pair<double,int> > ehdepth;
      spr::energyHCALCell(hotCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, depthHE, debug);
      for (unsigned int i=0; i<ehdepth.size(); ++i) {
        eHcalDepthHot[ehdepth[i].second-1] = ehdepth[i].first;
        HcalSubdetector subdet0 = (hbhe) ? ((ehdepth[i].second >= depthHE) ? HcalEndcap : HcalBarrel) : subdet;
        HcalDetId hcid0(subdet0,ieta,iphi,ehdepth[i].second);
        double actL = activeLength(DetId(hcid0));
        activeHotL[ehdepth[i].second-1] = actL;
        activeLengthHotTot += actL;
#ifdef EDM_ML_DEBUG
        if ((verbosity_%10) > 0)
          std::cout << hcid0 << " E " << ehdepth[i].first << " L " 
            << actL << std::endl;
#endif
      }
    }
#ifdef EDM_ML_DEBUG
    if ((verbosity_%10) > 0) {
      for (int k=0; k<depthMax_; ++k)
        std::cout << "Depth " << k << " E " << eHcalDepth[k] << ":" 
              << eHcalDepthHot[k] << std::endl;
    }
#endif     
    matchedId_.push_back(tmpmatch);                    
    ecal3x3Energy_.push_back(eEcal);
    ecalDetId_.push_back(isoCell.rawId());
    hcal1x1Energy_.push_back(eHcal);
    hcalDetId_.push_back(closestCell.rawId());
    for (int k=0; k<depthMax_; ++k) {
      hcalDepthEnergy_[k].push_back(eHcalDepth[k]);
      hcalDepthActiveLength_[k].push_back(activeL[k]);
      hcalDepthEnergyHot_[k].push_back(eHcalDepthHot[k]);
      hcalDepthActiveLengthHot_[k].push_back(activeHotL[k]);
    }
    hcalHot_.push_back(isHot);
    hcalActiveLengthHot_.push_back(activeLengthHotTot);
      }
    }
  }
  if (!hcalHot_.empty()) tree_->Fill();
}

void HcalHBHEMuonSimAnalyzer::beginJob() {

  tree_ = fs_->make<TTree>("TREE", "TREE");
  tree_->Branch("Run_No",           &runNumber_);
  tree_->Branch("Event_No",         &eventNumber_);
  tree_->Branch("LumiNumber",       &lumiNumber_);
  tree_->Branch("BXNumber",         &bxNumber_);
  tree_->Branch("pt_of_muon",       &ptGlob_);
  tree_->Branch("eta_of_muon",      &etaGlob_);
  tree_->Branch("phi_of_muon",      &phiGlob_);
  tree_->Branch("p_of_muon",        &pMuon_);
  tree_->Branch("matchedId",        &matchedId_);
  
  tree_->Branch("ecal_3x3",         &ecal3x3Energy_);
  tree_->Branch("ecal_detID",       &ecalDetId_);
  tree_->Branch("hcal_1x1",         &hcal1x1Energy_);
  tree_->Branch("hcal_detID",       &hcalDetId_);
  tree_->Branch("hcal_cellHot",     &hcalHot_);
  tree_->Branch("activeLength",     &hcalActiveLength_);
  tree_->Branch("activeLengthHot",  &hcalActiveLengthHot_);
  char name[100];
  for (int k=0; k<maxDepth_; ++k) {
    sprintf (name, "hcal_edepth%d", (k+1));
    tree_->Branch(name, &hcalDepthEnergy_[k]);
    sprintf (name, "hcal_activeL%d", (k+1));
    tree_->Branch(name,  &hcalDepthActiveLength_[k]);
    sprintf (name, "hcal_edepthHot%d", (k+1));
    tree_->Branch(name,  &hcalDepthEnergyHot_[k]);
    sprintf (name, "hcal_activeHotL%d", (k+1));
    tree_->Branch(name, &hcalDepthActiveLength_[k]);
  }
}

void HcalHBHEMuonSimAnalyzer::beginRun(edm::Run const& iRun, 
                       edm::EventSetup const& iSetup) {

  edm::ESHandle<HcalDDDRecConstants> pHRNDC;
  iSetup.get<HcalRecNumberingRecord>().get(pHRNDC);
  const HcalDDDRecConstants & hdc = (*pHRNDC);
  actHB_.clear();
  actHE_.clear();
  actHB_ = hdc.getThickActive(0);
  actHE_ = hdc.getThickActive(1);
}


void HcalHBHEMuonSimAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {

  edm::ParameterSetDescription desc;
  desc.add<std::string>("ModuleLabel","g4SimHits");
  desc.add<std::string>("EBCollection","EcalHitsEB");
  desc.add<std::string>("EECollection","EcalHitsEE");
  desc.add<std::string>("HCCollection","HcalHits");
  desc.addUntracked<int>("Verbosity",0);
  desc.addUntracked<int>("MaxDepth",4);
  desc.addUntracked<double>("EtaMax",3.0);
  desc.addUntracked<double>("TimeMinCutECAL",-500.);
  desc.addUntracked<double>("TimeMaxCutECAL",500.);
  desc.addUntracked<double>("TimeMinCutHCAL",-500.);
  desc.addUntracked<double>("TimeMaxCutHCAL",500.);
  descriptions.add("hcalHBHEMuonSim",desc);
}

void HcalHBHEMuonSimAnalyzer::clearVectors() {

  runNumber_   = -99999;
  eventNumber_ = -99999;
  lumiNumber_  = -99999;
  bxNumber_    = -99999;
    
  ptGlob_.clear();
  etaGlob_.clear(); 
  phiGlob_.clear(); 
  pMuon_.clear();
  matchedId_.clear();
  ecal3x3Energy_.clear();
  ecalDetId_.clear();
  hcal1x1Energy_.clear();
  hcalDetId_.clear();
  hcalHot_.clear();
  hcalActiveLength_.clear();
  hcalActiveLengthHot_.clear();
  for (int k=0; k<depthMax_; ++k) {
    hcalDepthEnergy_[k].clear();
    hcalDepthActiveLength_[k].clear();
    hcalDepthEnergyHot_[k].clear();
    hcalDepthActiveLengthHot_[k].clear();
  }
}

unsigned int HcalHBHEMuonSimAnalyzer::matchId(const HcalDetId& id1, 
                          const HcalDetId& id2) {
  
  HcalDetId kd1(id1.subdet(),id1.ieta(),id1.iphi(),1);
  HcalDetId kd2(id2.subdet(),id2.ieta(),id2.iphi(),1);
  unsigned int match = ((kd1 == kd2) ? 1 : 0);
  return match;
}

double HcalHBHEMuonSimAnalyzer::activeLength(const DetId& id_) {
  HcalDetId id(id_);
  int ieta = id.ietaAbs();
  int depth= id.depth();
  double lx(0);
  if (id.subdet() == HcalBarrel) {
    for (unsigned int i=0; i<actHB_.size(); ++i) {
      if (ieta == actHB_[i].ieta && depth == actHB_[i].depth) {
    lx = actHB_[i].thick;
    break;
      }
    }
  } else {
    for (unsigned int i=0; i<actHE_.size(); ++i) {
      if (ieta == actHE_[i].ieta && depth == actHE_[i].depth) {
    lx = actHE_[i].thick;
    break;
      }
    }
  }
  return lx;
}

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