AlCaHOCalibProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
// Dec 2015 : Added bool m_cosmic to choose cosmic or collision run through python file

//integrate the code with some 8_0_X or 7_6_X recent IB and run the
// following tests: 4.22, 8.0, 25.0, 140.53. You can always activate them using
//runTheMatrix.py -l 4.22

// 7th Nov 2015 :  tmpHOCalib.ecal03 = iso05.sumPt; // iso03.emEt+muonenr.em;
//                   tmpHOCalib.inslumi=lumiScale->begin()->pileup();
//
// April 2015 : Remove all digi part
//  Also look for HO geometry in CMSSW in parallel with stanalone one.
// Official one has problem in reco geometry, particularly tiles at the edge of wheel
// Remove all histogrammes except occupancy one
// Remove Trigger bits
// But addition of these variables, ilumi (analyser), inslumi (analyser), nprim

// Feb09 2009
// Move the initialisation of SteppingHelixPropagator from ::beginJob() to ::produce()
//
// Oct3 2008
// Difference in tag V00-02-45 with previous code

// 1. One new object on data format, which was realised in
//     CRUZET data analysis.
//2.  Remove all histogram and cout in the code
//3. An upgrade in code, which increases the acceptance of
//    muon near the edge (this also realised in CRUZET data).
// Difference in wrt V00-02-45
// 1. initialisation tmpHOCalib.htime = -1000;
// 2. By mistake HLT was commented out

// Package:    AlCaHOCalibProducer
// Class:      AlCaHOCalibProducer
//
//
// Original Author:  Gobinda Majumder
//         Created:  Fri Jul  6 17:17:21 CEST 2007
//
//

// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
#include "DataFormats/Luminosity/interface/LumiDetails.h"
#include "DataFormats/Math/interface/Error.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
#include "DataFormats/Scalers/interface/LumiScalers.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"

#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

// Necessary includes for identify severity of flagged problems in HO rechits
//#include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

#include "TH2F.h"

/* C++ Headers */
#include <string>

#include <iostream>
#include <fstream>
//
// class decleration
//

class AlCaHOCalibProducer : public edm::one::EDProducer<edm::one::SharedResources> {
public:
  explicit AlCaHOCalibProducer(const edm::ParameterSet&);
  ~AlCaHOCalibProducer() override = default;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
  typedef Basic3DVector<float> PositionType;
  typedef Basic3DVector<float> DirectionType;
  typedef Basic3DVector<float> RotationType;

private:
  void produce(edm::Event&, const edm::EventSetup&) override;
  void beginJob() override;
  void endJob() override;

  void fillHOStore(const reco::TrackRef& ncosm,
                   HOCalibVariables& tmpHOCalib,
                   std::unique_ptr<HOCalibVariableCollection>& hostore,
                   int Noccu_old,
                   int indx,
                   edm::Handle<reco::TrackCollection> cosmicmuon,
                   edm::View<reco::Muon>::const_iterator muon1,
                   const edm::Event& iEvent,
                   const CaloSubdetectorGeometry*,
                   const MagneticField&);

  void findHOEtaPhi(int iphsect, int& ietaho, int& iphiho);
  //  virtual void endRun(edm::Run const &, edm::EventSetup const &) override;
  // ----------member data ---------------------------

  float xhor0;  //x-position in ring 0
  float yhor0;  //y-position in ring 0
  float xhor1;  //x-position in ring 1
  float yhor1;  //y-position in ring 1
  int iring;    //Ring number -2,-1,0,1,2

  float localxhor0;  //local x-distance from edege in ring 0
  float localyhor0;  //local y-distance from edege in ring 0
  float localxhor1;  //local x-distance from edege in ring 1
  float localyhor1;  //local y-distance from edege in ring 1

  TH2F* ho_occupency[5];
  bool m_occupancy;
  bool m_cosmic;

  const int netabin = 16;
  const int nphimx = 72;
  const int netamx = 32;
  const int ncidmx = 5;
  const double rHOL0 = 382.0;
  const double rHOL1 = 407.0;

  edm::InputTag muonTags_;  // cosmicMuons (for cosmic run) or muons (for collision run)

  edm::EDGetTokenT<reco::TrackCollection> tok_muonsCosmic_;
  edm::EDGetTokenT<edm::View<reco::Muon> > tok_muons_;
  edm::EDGetTokenT<reco::VertexCollection> tok_vertex_;
  //  edm::EDGetTokenT<LumiDetails> tok_lumi_;
  edm::EDGetTokenT<LumiScalersCollection> tok_lumi_;

  edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
  edm::EDGetTokenT<HORecHitCollection> tok_ho_;
  edm::EDGetTokenT<CaloTowerCollection> tok_tower_;

  edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> tok_hcalChStatus_;
  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
  edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> tok_hcalSevLvlComputer_;
  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;

  bool m_hbinfo;
  int m_startTS;
  int m_endTS;
  double m_sigma;

  typedef math::Error<5>::type CovarianceMatrix;
  int Noccu;
  int nRuns;

  //  SteppingHelixPropagator* stepProp;
  FreeTrajectoryState getFreeTrajectoryState(const reco::Track& tk, const MagneticField* field, int itag, bool dir);

  unsigned int Ntp;  // # of HLT trigger paths (should be the same for all events!)
  std::map<std::string, bool> fired;

  //hcal severity ES
  const HcalChannelQuality* theHcalChStatus;
  const HcalSeverityLevelComputer* theHcalSevLvlComputer;
  int Nevents;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
AlCaHOCalibProducer::AlCaHOCalibProducer(const edm::ParameterSet& iConfig) {
  usesResource(TFileService::kSharedResource);
  //register your products

  m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
  m_sigma = iConfig.getUntrackedParameter<double>("sigma", 0.05);
  m_occupancy = iConfig.getUntrackedParameter<bool>("plotOccupancy", false);
  m_cosmic = iConfig.getUntrackedParameter<bool>("CosmicData", false);

  // keep InputTag muonTags_ since it is used below. - cowden
  muonTags_ = iConfig.getUntrackedParameter<edm::InputTag>("muons");
  tok_muonsCosmic_ = consumes<reco::TrackCollection>(muonTags_);
  tok_muons_ = consumes<edm::View<reco::Muon> >(muonTags_);
  tok_vertex_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexTags"));
  //  tok_lumi_ = consumes<LumiDetails ,edm::InLumi>(iConfig.getParameter<edm::InputTag>("lumiTags"));
  tok_lumi_ = consumes<LumiScalersCollection>(iConfig.getParameter<edm::InputTag>("lumiTags"));
  tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"));
  tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
  tok_tower_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("towerInput"));

  tok_hcalChStatus_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
  tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
  tok_hcalSevLvlComputer_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
  tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();

  produces<HOCalibVariableCollection>("HOCalibVariableCollection").setBranchAlias("HOCalibVariableCollection");

  if (m_occupancy) {
    edm::Service<TFileService> fs;

    char title[200];

    for (int ij = 0; ij < 5; ij++) {
      sprintf(title, "ho_occupency (>%i #sigma)", ij + 2);
      ho_occupency[ij] =
          fs->make<TH2F>(title, title, netamx + 1, -netamx - 0.5, netamx / 2 + 0.5, nphimx, 0.5, nphimx + 0.5);
    }
  }
}

//
// member functions
//

void AlCaHOCalibProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("hbheInput", edm::InputTag("hbhereco"));
  desc.addUntracked<bool>("hotime", false);
  desc.addUntracked<bool>("hbinfo", false);
  desc.addUntracked<double>("sigma", 1.0);
  desc.addUntracked<bool>("plotOccupancy", false);
  desc.addUntracked<bool>("CosmicData", false);
  desc.add<edm::InputTag>("hoInput", edm::InputTag("horeco"));
  desc.add<edm::InputTag>("towerInput", edm::InputTag("towerMaker"));
  desc.addUntracked<std::string>("RootFileName", "test.root");
  desc.addUntracked<double>("m_scale", 4.0);
  desc.addUntracked<bool>("debug", false);
  desc.addUntracked<edm::InputTag>("muons", edm::InputTag("muons"));
  desc.add<edm::InputTag>("vertexTags", edm::InputTag("offlinePrimaryVertices"));
  desc.add<edm::InputTag>("lumiTags", edm::InputTag("scalersRawToDigi"));
  descriptions.add("alcaHOCalibProducer", desc);
}

// ------------ method called to produce the data  ------------
void AlCaHOCalibProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  int irun = iEvent.id().run();
  //  int ilumi = iEvent.luminosityBlock();

  Nevents++;

  if (Nevents % 5000 == 1)
    edm::LogInfo("HOCalib") << "AlCaHOCalibProducer Processing event # " << Nevents << " " << Noccu << " " << irun
                            << " " << iEvent.id().event();

  theHcalChStatus = &iSetup.getData(tok_hcalChStatus_);

  auto hostore = std::make_unique<HOCalibVariableCollection>();

  edm::Handle<reco::TrackCollection> cosmicmuon;
  edm::Handle<edm::View<reco::Muon> > collisionmuon;

  bool muonOK(true);
  HOCalibVariables tmpHOCalib;
  tmpHOCalib.nprim = -1;
  tmpHOCalib.inslumi = -1.;

  if (m_cosmic) {
    cosmicmuon = iEvent.getHandle(tok_muonsCosmic_);
    muonOK = (cosmicmuon.isValid() && !cosmicmuon->empty());
  } else {
    collisionmuon = iEvent.getHandle(tok_muons_);
    muonOK = (collisionmuon.isValid() && !collisionmuon->empty());

    if (iEvent.isRealData()) {
      auto const& primaryVertices = iEvent.getHandle(tok_vertex_);
      if (primaryVertices.isValid()) {
        tmpHOCalib.nprim = primaryVertices->size();
      }

      tmpHOCalib.inslumi = 0.;

      auto const& lumiScale = iEvent.getHandle(tok_lumi_);

      if (lumiScale.isValid()) {
        if (lumiScale->empty()) {
          edm::LogError("HOCalib") << "lumiScale collection is empty";
        } else {
          tmpHOCalib.inslumi = lumiScale->begin()->pileup();
        }
      }
    }
  }

  if (muonOK) {
    int Noccu_old = Noccu;
    edm::View<reco::Muon>::const_iterator muon1;

    theHcalSevLvlComputer = &iSetup.getData(tok_hcalSevLvlComputer_);

    MagneticField const& magField = iSetup.getData(tok_magField_);

    const CaloGeometry& geo = iSetup.getData(tok_geom_);
    const CaloSubdetectorGeometry* gHO = geo.getSubdetectorGeometry(DetId::Hcal, HcalOuter);

    if (m_cosmic) {
      int indx(0);
      for (reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin(); ncosm != cosmicmuon->end();
           ++ncosm, ++indx) {
        if ((*ncosm).ndof() < 15)
          continue;
        if ((*ncosm).normalizedChi2() > 30.0)
          continue;
        reco::TrackRef tRef = reco::TrackRef(cosmicmuon, indx);
        fillHOStore(tRef, tmpHOCalib, hostore, Noccu_old, indx, cosmicmuon, muon1, iEvent, gHO, magField);
      }
    } else {
      for (muon1 = collisionmuon->begin(); muon1 < collisionmuon->end(); muon1++) {
        if ((!muon1->isGlobalMuon()) || (!muon1->isTrackerMuon()))
          continue;
        reco::TrackRef ncosm = muon1->innerTrack();
        fillHOStore(ncosm, tmpHOCalib, hostore, Noccu_old, 0, cosmicmuon, muon1, iEvent, gHO, magField);
      }
    }
  }

  iEvent.put(std::move(hostore), "HOCalibVariableCollection");
}

// ------------ method called once each job just before starting event loop  ------------
void AlCaHOCalibProducer::beginJob() {
  Nevents = 0;
  nRuns = 0;
  Noccu = 0;
}

// ------------ method called once each job just after ending the event loop  ------------
void AlCaHOCalibProducer::endJob() {
  if (m_occupancy) {
    for (int ij = 0; ij < 5; ij++) {
      ho_occupency[ij]->Scale(1. / std::max(1, Noccu));
    }
  }
  edm::LogInfo("HOCalib") << " AlCaHOCalibProducer processed event " << Nevents;
}

void AlCaHOCalibProducer::fillHOStore(const reco::TrackRef& ncosm,
                                      HOCalibVariables& tmpHOCalib,
                                      std::unique_ptr<HOCalibVariableCollection>& hostore,
                                      int Noccu_old,
                                      int indx,
                                      edm::Handle<reco::TrackCollection> cosmicmuon,
                                      edm::View<reco::Muon>::const_iterator muon1,
                                      const edm::Event& iEvent,
                                      const CaloSubdetectorGeometry* gHO,
                                      const MagneticField& magField) {
  // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined

  int charge = ncosm->charge();

  double innerr = (*ncosm).innerPosition().Perp2();
  double outerr = (*ncosm).outerPosition().Perp2();
  int iiner = (innerr < outerr) ? 1 : 0;

  //---------------------------------------------------
  //             in_to_out  Dir         in_to_out  Dir
  //   StandAlone ^         ^     Cosmic    ^    |
  //              |         |               |    v
  //---------------------------------------------------Y=0
  //   StandAlone |         |     Cosmic    ^    |
  //              v         v               |    v
  //----------------------------------------------------

  double posx, posy, posz;
  double momx, momy, momz;

  if (iiner == 1) {
    posx = (*ncosm).innerPosition().X();
    posy = (*ncosm).innerPosition().Y();
    posz = (*ncosm).innerPosition().Z();

    momx = (*ncosm).innerMomentum().X();
    momy = (*ncosm).innerMomentum().Y();
    momz = (*ncosm).innerMomentum().Z();

  } else {
    posx = (*ncosm).outerPosition().X();
    posy = (*ncosm).outerPosition().Y();
    posz = (*ncosm).outerPosition().Z();

    momx = (*ncosm).outerMomentum().X();
    momy = (*ncosm).outerMomentum().Y();
    momz = (*ncosm).outerMomentum().Z();
  }

  PositionType trkpos(posx, posy, posz);

  CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
  CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);

  bool samedir = (tmpmuon3v.dot(tmpmuondir) > 0) ? true : false;
  for (int ij = 0; ij < 3; ij++) {
    tmpHOCalib.caloen[ij] = 0.0;
  }
  int inearbymuon = 0;
  localxhor0 = localyhor0 = 20000;  //GM for 22OCT07 data

  if (m_cosmic) {
    int ind(0);
    for (reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin(); ncosmcor != cosmicmuon->end();
         ++ncosmcor, ++ind) {
      if (indx == ind)
        continue;
      CLHEP::Hep3Vector tmpmuon3vcor;
      CLHEP::Hep3Vector tmpmom3v;
      if (iiner == 1) {
        tmpmuon3vcor = CLHEP::Hep3Vector(
            (*ncosmcor).innerPosition().X(), (*ncosmcor).innerPosition().Y(), (*ncosmcor).innerPosition().Z());
        tmpmom3v = CLHEP::Hep3Vector(
            (*ncosmcor).innerMomentum().X(), (*ncosmcor).innerMomentum().Y(), (*ncosmcor).innerMomentum().Z());
      } else {
        tmpmuon3vcor = CLHEP::Hep3Vector(
            (*ncosmcor).outerPosition().X(), (*ncosmcor).outerPosition().Y(), (*ncosmcor).outerPosition().Z());
        tmpmom3v = CLHEP::Hep3Vector(
            (*ncosmcor).outerMomentum().X(), (*ncosmcor).outerMomentum().Y(), (*ncosmcor).outerMomentum().Z());
      }

      if (tmpmom3v.mag() < 0.2 || (*ncosmcor).ndof() < 5)
        continue;

      double angle = tmpmuon3v.angle(tmpmuon3vcor);
      if (angle < 7.5 * CLHEP::deg) {
        inearbymuon = 1;
      }  //  break;}

      //    if (muonTagsi_.label() =="cosmicMuons") {
      if (angle < 7.5 * CLHEP::deg) {
        tmpHOCalib.caloen[0] += 1.;
      }
      if (angle < 15.0 * CLHEP::deg) {
        tmpHOCalib.caloen[1] += 1.;
      }
      if (angle < 35.0 * CLHEP::deg) {
        tmpHOCalib.caloen[2] += 1.;
      }
    }
  } else {
    //            if (muonTags_.label() =="muons") {
    auto const& calotower = iEvent.getHandle(tok_tower_);

    for (CaloTowerCollection::const_iterator calt = calotower->begin(); calt != calotower->end(); calt++) {
      //CMSSW_2_1_x const math::XYZVector towermom = (*calt).momentum();
      double ith = (*calt).momentum().theta();
      double iph = (*calt).momentum().phi();

      CLHEP::Hep3Vector calo3v(sin(ith) * cos(iph), sin(ith) * sin(iph), cos(ith));

      double angle = tmpmuon3v.angle(calo3v);

      if (angle < 7.5 * CLHEP::deg) {
        tmpHOCalib.caloen[0] += calt->emEnergy() + calt->hadEnergy();
      }
      if (angle < 15 * CLHEP::deg) {
        tmpHOCalib.caloen[1] += calt->emEnergy() + calt->hadEnergy();
      }
      if (angle < 35 * CLHEP::deg) {
        tmpHOCalib.caloen[2] += calt->emEnergy() + calt->hadEnergy();
      }
    }
  }
  if ((m_cosmic) || (tmpHOCalib.caloen[0] <= 10.0)) {
    GlobalPoint glbpt(posx, posy, posz);

    double mom = sqrt(momx * momx + momy * momy + momz * momz);

    momx /= mom;
    momy /= mom;
    momz /= mom;

    DirectionType trkdir(momx, momy, momz);

    tmpHOCalib.trkdr = (*ncosm).d0();
    tmpHOCalib.trkdz = (*ncosm).dz();
    tmpHOCalib.nmuon = (m_cosmic) ? cosmicmuon->size() : 1;
    tmpHOCalib.trkvx = glbpt.x();
    tmpHOCalib.trkvy = glbpt.y();
    tmpHOCalib.trkvz = glbpt.z();
    tmpHOCalib.trkmm = mom * charge;
    tmpHOCalib.trkth = trkdir.theta();
    tmpHOCalib.trkph = trkdir.phi();
    tmpHOCalib.isect2 = -2;
    tmpHOCalib.isect = -2;
    tmpHOCalib.hodx = -100;
    tmpHOCalib.hody = -100;
    tmpHOCalib.hoang = -2.0;
    tmpHOCalib.momatho = -2;
    tmpHOCalib.ndof = (inearbymuon == 0) ? (int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
    tmpHOCalib.chisq = (*ncosm).normalizedChi2();  // max(1.,tmpHOCalib.ndof);
    if (!m_cosmic) {
      reco::MuonEnergy muonenr = muon1->calEnergy();
      reco::MuonIsolation iso03 = muon1->isolationR03();
      reco::MuonIsolation iso05 = muon1->isolationR05();

      tmpHOCalib.tkpt03 = iso03.sumPt;
      tmpHOCalib.ecal03 = iso05.sumPt;  // iso03.emEt+muonenr.em;
      tmpHOCalib.hcal03 = iso03.hadEt + muonenr.had;
    }
    tmpHOCalib.therr = 0.;
    tmpHOCalib.pherr = 0.;
    if (iiner == 1) {
      reco::TrackBase::CovarianceMatrix innercov = (*ncosm).innerStateCovariance();
      tmpHOCalib.therr = innercov(1, 1);  //thetaError();
      tmpHOCalib.pherr = innercov(2, 2);  //phi0Error();
    } else {
      reco::TrackBase::CovarianceMatrix outercov = (*ncosm).outerStateCovariance();
      tmpHOCalib.therr = outercov(1, 1);  //thetaError();
      tmpHOCalib.pherr = outercov(2, 2);  //phi0Error();
    }

    SteppingHelixPropagator myHelix(&magField, anyDirection);
    myHelix.setMaterialMode(false);
    myHelix.applyRadX0Correction(true);
    double phiho = trkpos.phi();
    if (phiho < 0)
      phiho += CLHEP::twopi;

    int iphisect_dt = int(6 * (phiho + 10.0 * CLHEP::deg) / CLHEP::pi);  //for u 18/12/06
    if (iphisect_dt >= 12)
      iphisect_dt = 0;

    int iphisect = -1;
    bool ipath = false;
    for (int kl = 0; kl <= 2; kl++) {
      int iphisecttmp = (kl < 2) ? iphisect_dt + kl : iphisect_dt - 1;
      if (iphisecttmp < 0)
        iphisecttmp = 11;
      if (iphisecttmp >= 12)
        iphisecttmp = 0;

      double phipos = iphisecttmp * CLHEP::pi / 6.;
      double phirot = phipos;

      GlobalVector xLocal(-sin(phirot), cos(phirot), 0.);
      GlobalVector yLocal(0., 0., 1.);
      GlobalVector zLocal = xLocal.cross(yLocal).unit();
      //    GlobalVector zLocal(cos(phirot), sin(phirot), 0.0);

      FreeTrajectoryState freetrajectorystate_ = getFreeTrajectoryState(*ncosm, &(magField), iiner, samedir);

      Surface::RotationType rot(xLocal, yLocal, zLocal);

      for (int ik = 1; ik >= 0; ik--) {  //propagate track in two HO layers

        double radial = rHOL1;
        if (ik == 0)
          radial = rHOL0;

        Surface::PositionType pos(radial * cos(phipos), radial * sin(phipos), 0.);
        PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);

        auto aPlane2 = new Plane(pos, rot);

        SteppingHelixStateInfo steppingHelixstateinfo_;
        myHelix.propagate(SteppingHelixStateInfo(freetrajectorystate_), (*aPlane2), steppingHelixstateinfo_);

        if (steppingHelixstateinfo_.isValid()) {
          GlobalPoint hotrkpos2xx(steppingHelixstateinfo_.position().x(),
                                  steppingHelixstateinfo_.position().y(),
                                  steppingHelixstateinfo_.position().z());

          if (ik == 1) {
            HcalDetId ClosestCell = (HcalDetId)gHO->getClosestCell(hotrkpos2xx);
            int ixeta = ClosestCell.ieta();
            int ixphi = ClosestCell.iphi();
            tmpHOCalib.isect2 = 100 * std::abs(ixeta + 50) + std::abs(ixphi);
          }

          GlobalVector hotrkpos2(steppingHelixstateinfo_.position().x(),
                                 steppingHelixstateinfo_.position().y(),
                                 steppingHelixstateinfo_.position().z());
          CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.momentum().x(),
                                      steppingHelixstateinfo_.momentum().y(),
                                      steppingHelixstateinfo_.momentum().z());

          LocalVector lclvt0 = (*aPlane).toLocal(hotrkpos2);

          double xx = lclvt0.x();
          double yy = lclvt0.y();

          if (ik == 1) {
            if ((std::abs(yy) < 130 && xx > -64.7 && xx < 138.2)                              //Ring-0
                || (std::abs(yy) > 130 && std::abs(yy) < 700 && xx > -76.3 && xx < 140.5)) {  //Ring +-1,2
              ipath = true;  //Only look for tracks which as hits in layer 1
              iphisect = iphisecttmp;
            }
          }

          if (iphisect != iphisecttmp)
            continue;  //Look for ring-0 only when ring1 is accepted for that sector

          switch (ik) {
            case 0:
              xhor0 = xx;  //lclvt0.x();
              yhor0 = yy;  //lclvt0.y();
              break;
            case 1:
              xhor1 = xx;  //lclvt0.x();
              yhor1 = yy;  //lclvt0.y();
              tmpHOCalib.momatho = hotrkdir2.mag();
              tmpHOCalib.hoang = CLHEP::Hep3Vector(zLocal.x(), zLocal.y(), zLocal.z()).dot(hotrkdir2.unit());
              break;
            default:
              break;
          }
        } else {
          break;
        }
      }
      if (ipath)
        break;
    }
    if (ipath) {  //If muon crossed HO laeyrs

      int ietaho = 50;
      int iphiho = -1;

      for (int ij = 0; ij < 9; ij++) {
        tmpHOCalib.hosig[ij] = -100.0;
      }
      for (int ij = 0; ij < 18; ij++) {
        tmpHOCalib.hocorsig[ij] = -100.0;
      }
      for (int ij = 0; ij < 9; ij++) {
        tmpHOCalib.hbhesig[ij] = -100.0;
      }
      tmpHOCalib.hocro = -100;
      tmpHOCalib.htime = -1000;

      int isect = 0;

      findHOEtaPhi(iphisect, ietaho, iphiho);

      if (ietaho != 0 && iphiho != 0 && std::abs(iring) <= 2) {  //Muon passed through a tower
        isect = 100 * std::abs(ietaho + 50) + std::abs(iphiho);
        if (std::abs(ietaho) >= netabin || iphiho < 0)
          isect *= -1;  //Not extrapolated to any tower
        if (std::abs(ietaho) >= netabin)
          isect -= 1000000;  //not matched with eta
        if (iphiho < 0)
          isect -= 2000000;  //not matched with phi
        tmpHOCalib.isect = isect;

        tmpHOCalib.hodx = localxhor1;
        tmpHOCalib.hody = localyhor1;

        if (iring == 0) {
          tmpHOCalib.hocorsig[8] = localxhor0;
          tmpHOCalib.hocorsig[9] = localyhor0;
        }

        int etamn = -4;
        int etamx = 4;
        if (iring == 1) {
          etamn = 5;
          etamx = 10;
        }
        if (iring == 2) {
          etamn = 11;
          etamx = 16;
        }
        if (iring == -1) {
          etamn = -10;
          etamx = -5;
        }
        if (iring == -2) {
          etamn = -16;
          etamx = -11;
        }

        int phimn = 1;
        int phimx = 2;
        if (iring == 0) {
          phimx = 2 * int((iphiho + 1) / 2.);
          phimn = phimx - 1;
        } else {
          phimn = 3 * int((iphiho + 1) / 3.) - 1;
          phimx = phimn + 2;
        }

        if (phimn < 1)
          phimn += nphimx;
        if (phimx > 72)
          phimx -= nphimx;

        if (m_hbinfo) {
          for (int ij = 0; ij < 9; ij++) {
            tmpHOCalib.hbhesig[ij] = -100.0;
          }

          auto const& hbheht = iEvent.getHandle(tok_hbhe_);  // iEvent.getByType(hbheht);
          if (!(*hbheht).empty()) {
            if ((*hbheht).empty())
              throw(int)(*hbheht).size();

            for (HBHERecHitCollection::const_iterator jk = (*hbheht).begin(); jk != (*hbheht).end(); jk++) {
              HcalDetId id = (*jk).id();
              int tmpeta = id.ieta();
              int tmpphi = id.iphi();

              int deta = tmpeta - ietaho;
              if (tmpeta < 0 && ietaho > 0)
                deta += 1;
              if (tmpeta > 0 && ietaho < 0)
                deta -= 1;

              //        if (tmpeta==-1 && ietaho== 1) deta = -1;
              //        if (tmpeta== 1 && ietaho==-1) deta =  1;

              int dphi = tmpphi - iphiho;
              if (dphi > nphimx / 2) {
                dphi -= nphimx;
              }
              if (dphi < -nphimx / 2) {
                dphi += nphimx;
              }

              //        if (phimn >phimx) {
              //          if (dphi==71) dphi=-1;
              //          if (dphi==-71) dphi=1;
              //        }

              if (m_occupancy) {
                float signal = (*jk).energy();
                //      int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
                if (signal > -100 && Noccu == Noccu_old) {
                  for (int ij = 0; ij < 5; ij++) {
                    if (signal > (ij + 2) * m_sigma) {
                      ho_occupency[ij]->Fill(tmpeta, tmpphi);
                    }
                  }
                }
              }

              int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;  //NEED correction in full CMS detector
              if (ipass2 == 0)
                continue;

              float signal = (*jk).energy();

              if (3 * (deta + 1) + dphi + 1 < 9)
                tmpHOCalib.hbhesig[3 * (deta + 1) + dphi + 1] = signal;
            }
          }
        }  //m_hbinfo #endif

        auto const& hoht = iEvent.getHandle(tok_ho_);

        if (!(*hoht).empty()) {
          for (HORecHitCollection::const_iterator jk = (*hoht).begin(); jk != (*hoht).end(); jk++) {
            HcalDetId id = (*jk).id();
            int tmpeta = id.ieta();
            int tmpphi = id.iphi();

            int ipass1 = 0;
            if (tmpeta >= etamn && tmpeta <= etamx) {
              if (phimn < phimx) {
                ipass1 = (tmpphi >= phimn && tmpphi <= phimx) ? 1 : 0;
              } else {
                ipass1 = (tmpphi == 71 || tmpphi == 72 || tmpphi == 1) ? 1 : 0;
              }
            }

            int deta = tmpeta - ietaho;
            int dphi = tmpphi - iphiho;

            if (tmpeta < 0 && ietaho > 0)
              deta += 1;
            if (tmpeta > 0 && ietaho < 0)
              deta -= 1;
            //        if (tmpeta==-1 && ietaho== 1) deta = -1;
            //        if (tmpeta== 1 && ietaho==-1) deta =  1;

            if (dphi > nphimx / 2) {
              dphi -= nphimx;
            }
            if (dphi < -nphimx / 2) {
              dphi += nphimx;
            }
            //        if (phimn>phimx) {
            //      if (dphi==71) dphi=-1;
            //      if (dphi==-71) dphi=1;
            //        }

            float signal = (*jk).energy();

            int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;

            if (ipass1 == 0 && ipass2 == 0)
              continue;

            if (ipass1 == 1) {
              int tmpdph = tmpphi - phimn;
              if (tmpdph < 0)
                tmpdph = 2;  //only case of iphi==1, where phimn=71

              int ilog = 2 * (tmpeta - etamn) + tmpdph;
              if (iring != 0) {
                if (iring > 0) {
                  ilog = 3 * (tmpeta - etamn) + tmpdph;  //Again CMS correction
                } else {
                  ilog = 3 * (etamx - tmpeta) + tmpdph;  //Again CMS correction
                }
              }
              if (ilog > -1 && ilog < 18) {
                tmpHOCalib.hocorsig[ilog] = signal;
              }
            }

            if (ipass2 == 1) {
              if (3 * (deta + 1) + dphi + 1 < 9) {
                tmpHOCalib.hosig[3 * (deta + 1) + dphi + 1] = signal;  //Again CMS azimuthal near phi 1&72
              }
            }

            if (deta == 0 && dphi == 0) {
              tmpHOCalib.htime = (*jk).time();
              tmpHOCalib.hoflag = (*jk).flags();

              // Get Channel Quality information for the given detID
              unsigned theStatusValue = theHcalChStatus->getValues(id)->getValue();
              // Now get severity of problems for the given detID, based on the rechit flag word and the channel quality status value
              int hitSeverity = theHcalSevLvlComputer->getSeverityLevel(id, (*jk).flags(), theStatusValue);
              tmpHOCalib.hoflag = hitSeverity;
              int crphi = tmpphi + 6;
              if (crphi > 72)
                crphi -= 72;

              for (HORecHitCollection::const_iterator jcr = (*hoht).begin(); jcr != (*hoht).end(); jcr++) {
                const HORecHit reccr = (const HORecHit)(*jcr);
                HcalDetId idcr = reccr.id();
                int etacr = idcr.ieta();
                int phicr = idcr.iphi();
                if (tmpeta == etacr && crphi == phicr) {
                  tmpHOCalib.hocro = reccr.energy();
                }
              }
            }
          }
        }
      }

      //GMA   Npass++;
      if (Noccu == Noccu_old)
        Noccu++;
      hostore->push_back(tmpHOCalib);
    }  // if (ipath)
  }    // Cut on calo energy
}

void AlCaHOCalibProducer::findHOEtaPhi(int iphisect, int& ietaho, int& iphiho) {
  //18/12/06 : use only position, not angle phi

  const double etalow[16] = {0.025,
                             35.195,
                             70.625,
                             106.595,
                             141.565,
                             180.765,
                             220.235,
                             261.385,
                             304.525,
                             349.975,
                             410.025,
                             452.085,
                             506.645,
                             565.025,
                             627.725,
                             660.25};
  const double etahgh[16] = {35.145,
                             70.575,
                             106.545,
                             125.505,
                             180.715,
                             220.185,
                             261.335,
                             304.475,
                             349.925,
                             392.575,
                             452.035,
                             506.595,
                             564.975,
                             627.675,
                             661.075,
                             700.25};

  const double philow[6] = {-76.27, -35.11, 0.35, 35.81, 71.77, 108.93};  //Ring+/-1 & 2
  const double phihgh[6] = {-35.81, -0.35, 35.11, 71.07, 108.23, 140.49};

  const double philow00[6] = {-60.27, -32.91, 0.35, 33.61, 67.37, 102.23};  //Ring0 L0
  const double phihgh00[6] = {-33.61, -0.35, 32.91, 66.67, 101.53, 129.49};

  const double philow01[6] = {-64.67, -34.91, 0.35, 35.61, 71.37, 108.33};  //Ring0 L1
  const double phihgh01[6] = {-35.61, -0.35, 34.91, 70.67, 107.63, 138.19};

  iring = -10;

  double tmpdy = std::abs(yhor1);
  for (int ij = 0; ij < netabin; ij++) {
    if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
      ietaho = ij + 1;
      float tmp1 = fabs(tmpdy - etalow[ij]);
      float tmp2 = fabs(tmpdy - etahgh[ij]);

      localyhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
      if (yhor1 < 0)
        localyhor1 *= -1.;

      if (ij < 4)
        iring = 0;
      if (ij >= 4 && ij < 10)
        iring = 1;
      if (ij >= 10 && ij < netabin)
        iring = 2;
      break;
    }
  }

  int tmpphi = 0;
  int tmpphi0 = 0;

  if (ietaho > 4) {  //Ring 1 and 2
    for (int ij = 0; ij < 6; ij++) {
      if (xhor1 > philow[ij] && xhor1 < phihgh[ij]) {
        tmpphi = ij + 1;
        float tmp1 = fabs(xhor1 - philow[ij]);
        float tmp2 = fabs(xhor1 - phihgh[ij]);
        localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
        break;
      }
    }
  } else {  //Ring 0
    for (int ij = 0; ij < 6; ij++) {
      if (xhor1 > philow01[ij] && xhor1 < phihgh01[ij]) {
        tmpphi = ij + 1;
        float tmp1 = fabs(xhor1 - philow01[ij]);
        float tmp2 = fabs(xhor1 - phihgh01[ij]);
        localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
        break;
      }
    }

    for (int ij = 0; ij < 6; ij++) {
      if (xhor0 > philow00[ij] && xhor0 < phihgh00[ij]) {
        tmpphi0 = ij + 1;
        float tmp1 = fabs(xhor0 - philow00[ij]);
        float tmp2 = fabs(xhor0 - phihgh00[ij]);
        localxhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
        if (tmpphi != tmpphi0)
          localxhor0 += 10000.;
        break;
      }
    }

    double tmpdy = std::abs(yhor0);
    for (int ij = 0; ij < 4; ij++) {
      if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
        float tmp1 = fabs(tmpdy - etalow[ij]);
        float tmp2 = fabs(tmpdy - etahgh[ij]);
        localyhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
        if (yhor0 < 0)
          localyhor0 *= -1.;
        if (ij + 1 != ietaho)
          localyhor0 += 10000.;
        break;
      }
    }
  }

  if (tmpphi != 0) {
    iphiho = 6 * iphisect - 2 + tmpphi;
    if (iphiho <= 0)
      iphiho += nphimx;
    if (iphiho > nphimx)
      iphiho -= nphimx;
  }

  //  isect2 = 15*iring+iphisect+1;

  if (yhor1 < 0) {
    if (std::abs(ietaho) > netabin) {  //Initialised with 50
      ietaho += 1;
    } else {
      ietaho *= -1;
    }
    //    isect2 *=-1;
    iring *= -1;
  }
}

FreeTrajectoryState AlCaHOCalibProducer::getFreeTrajectoryState(const reco::Track& tk,
                                                                const MagneticField* field,
                                                                int iiner,
                                                                bool dir) {
  if (iiner == 0) {
    GlobalPoint gpos(tk.outerX(), tk.outerY(), tk.outerZ());
    GlobalVector gmom(tk.outerPx(), tk.outerPy(), tk.outerPz());
    if (dir)
      gmom *= -1.;
    GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
    CurvilinearTrajectoryError err(tk.extra()->outerStateCovariance());
    return FreeTrajectoryState(par, err);
  } else {
    GlobalPoint gpos(tk.innerPosition().X(), tk.innerPosition().Y(), tk.innerPosition().Z());
    GlobalVector gmom(tk.innerMomentum().X(), tk.innerMomentum().Y(), tk.innerMomentum().Z());
    if (dir)
      gmom *= -1.;
    GlobalTrajectoryParameters par(gpos, -gmom, tk.charge(), field);
    CurvilinearTrajectoryError err(tk.extra()->innerStateCovariance());
    return FreeTrajectoryState(par, err);
  }
}

#include "FWCore/Framework/interface/MakerMacros.h"

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