CastorClusterProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    Castor
// Class:      CastorClusterProducer
//
//
// Original Author:  Hans Van Haevermaet, Benoit Roland
//         Created:  Wed Jul  9 14:00:40 CEST 2008
//
//

// system include
#include <memory>
#include <vector>
#include <iostream>
#include <TMath.h>
#include <TRandom3.h>

// user include
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/Math/interface/Point3D.h"

// Castor Object include
#include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
#include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
#include "DataFormats/CastorReco/interface/CastorTower.h"
#include "DataFormats/CastorReco/interface/CastorCluster.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/JetReco/interface/BasicJet.h"
#include "DataFormats/JetReco/interface/BasicJetCollection.h"
#include "DataFormats/JetReco/interface/Jet.h"

//
// class decleration
//

class CastorClusterProducer : public edm::global::EDProducer<> {
public:
  explicit CastorClusterProducer(const edm::ParameterSet&);
  ~CastorClusterProducer() override;

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

  // ----------member data ---------------------------
  typedef math::XYZPointD Point;
  typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
  typedef ROOT::Math::RhoZPhiPoint CellPoint;
  typedef std::vector<reco::CastorTower> CastorTowerCollection;
  typedef std::vector<reco::CastorCluster> CastorClusterCollection;
  std::string input_, basicjets_;
  edm::EDGetTokenT<CastorTowerCollection> tok_input_;
  edm::EDGetTokenT<reco::BasicJetCollection> tok_jets_;
  edm::EDGetTokenT<CastorTowerCollection> tok_tower_;
  bool clusteralgo_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructor and destructor
//

CastorClusterProducer::CastorClusterProducer(const edm::ParameterSet& iConfig)
    : input_(iConfig.getUntrackedParameter<std::string>("inputtowers", "")),
      basicjets_(iConfig.getUntrackedParameter<std::string>("basicjets", "")),
      clusteralgo_(iConfig.getUntrackedParameter<bool>("ClusterAlgo", false)) {
  // register for data access
  tok_input_ = consumes<CastorTowerCollection>(edm::InputTag(input_));
  tok_jets_ = consumes<reco::BasicJetCollection>(edm::InputTag(basicjets_));
  tok_tower_ = consumes<CastorTowerCollection>(edm::InputTag("CastorTowerReco"));
  // register your products
  produces<CastorClusterCollection>();
  // now do what ever other initialization is needed
}

CastorClusterProducer::~CastorClusterProducer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//
// member functions
//

// ------------ method called to produce the data  ------------
void CastorClusterProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  using namespace edm;
  using namespace reco;
  using namespace TMath;

  LogDebug("CastorClusterProducer") << "3. entering CastorClusterProducer";

  if (!input_.empty()) {
    // Produce CastorClusters from CastorTowers

    edm::Handle<CastorTowerCollection> InputTowers;
    iEvent.getByToken(tok_input_, InputTowers);

    auto OutputClustersfromClusterAlgo = std::make_unique<CastorClusterCollection>();

    // get and check input size
    int nTowers = InputTowers->size();

    if (nTowers == 0)
      LogDebug("CastorClusterProducer") << "Warning: You are trying to run the Cluster algorithm with 0 input towers.";

    CastorTowerRefVector posInputTowers, negInputTowers;

    for (size_t i = 0; i < InputTowers->size(); ++i) {
      reco::CastorTowerRef tower_p = reco::CastorTowerRef(InputTowers, i);
      if (tower_p->eta() > 0.)
        posInputTowers.push_back(tower_p);
      if (tower_p->eta() < 0.)
        negInputTowers.push_back(tower_p);
    }

    // build cluster from ClusterAlgo
    if (clusteralgo_ == true) {
      // code
      iEvent.put(std::move(OutputClustersfromClusterAlgo));
    }
  }

  if (!basicjets_.empty()) {
    Handle<BasicJetCollection> bjCollection;
    iEvent.getByToken(tok_jets_, bjCollection);

    Handle<CastorTowerCollection> ctCollection;
    iEvent.getByToken(tok_tower_, ctCollection);

    auto OutputClustersfromBasicJets = std::make_unique<CastorClusterCollection>();

    if (bjCollection->empty())
      LogDebug("CastorClusterProducer")
          << "Warning: You are trying to run the Cluster algorithm with 0 input basicjets.";

    for (unsigned i = 0; i < bjCollection->size(); i++) {
      const BasicJet* bj = &(*bjCollection)[i];

      double energy = bj->energy();
      TowerPoint temp(88.5, bj->eta(), bj->phi());
      Point position(temp);
      double emEnergy = 0.;
      double hadEnergy = 0.;
      double width = 0.;
      double depth = 0.;
      double fhot = 0.;
      double sigmaz = 0.;
      CastorTowerRefVector usedTowers;
      double zmean = 0.;
      double z2mean = 0.;

      std::vector<CandidatePtr> ccp = bj->getJetConstituents();
      std::vector<CandidatePtr>::const_iterator itParticle;
      for (itParticle = ccp.begin(); itParticle != ccp.end(); ++itParticle) {
        const CastorTower* castorcand = dynamic_cast<const CastorTower*>(itParticle->get());
        //cout << " castortowercandidate reference energy = " << castorcand->castorTower()->energy() << endl;
        //cout << " castortowercandidate reference eta = " << castorcand->castorTower()->eta() << endl;
        //cout << " castortowercandidate reference phi = " << castorcand->castorTower()->phi() << endl;
        //cout << " castortowercandidate reference depth = " << castorcand->castorTower()->depth() << endl;

        //CastorTowerCollection *ctc = new CastorTowerCollection();
        //ctc->push_back(*castorcand);
        //CastorTowerRef towerref = CastorTowerRef(ctc,0);

        size_t thisone = 0;
        for (size_t l = 0; l < ctCollection->size(); l++) {
          const CastorTower ct = (*ctCollection)[l];
          if (std::abs(ct.phi() - castorcand->phi()) < 0.0001) {
            thisone = l;
          }
        }

        CastorTowerRef towerref(ctCollection, thisone);
        usedTowers.push_back(towerref);
        emEnergy += castorcand->emEnergy();
        hadEnergy += castorcand->hadEnergy();
        depth += castorcand->depth() * castorcand->energy();
        width += pow(phiangle(castorcand->phi() - bj->phi()), 2) * castorcand->energy();
        fhot += castorcand->fhot() * castorcand->energy();

        // loop over rechits
        for (edm::RefVector<edm::SortedCollection<CastorRecHit> >::iterator it = castorcand->rechitsBegin();
             it != castorcand->rechitsEnd();
             it++) {
          edm::Ref<edm::SortedCollection<CastorRecHit> > rechit_p = *it;
          double Erechit = rechit_p->energy();
          HcalCastorDetId id = rechit_p->id();
          int module = id.module();
          double zrechit = 0;
          if (module < 3)
            zrechit = -14390 - 24.75 - 49.5 * (module - 1);
          if (module > 2)
            zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
          zmean += Erechit * zrechit;
          z2mean += Erechit * zrechit * zrechit;
        }  // end loop over rechits
      }
      //cout << "" << endl;

      depth = depth / energy;
      width = sqrt(width / energy);
      fhot = fhot / energy;

      zmean = zmean / energy;
      z2mean = z2mean / energy;
      double sigmaz2 = z2mean - zmean * zmean;
      if (sigmaz2 > 0)
        sigmaz = sqrt(sigmaz2);

      CastorCluster cc(
          energy, position, emEnergy, hadEnergy, emEnergy / energy, width, depth, fhot, sigmaz, usedTowers);
      OutputClustersfromBasicJets->push_back(cc);
    }

    iEvent.put(std::move(OutputClustersfromBasicJets));
  }
}

// help function to calculate phi within [-pi,+pi]
double CastorClusterProducer::phiangle(double testphi) const {
  double phi = testphi;
  while (phi > M_PI)
    phi -= (2 * M_PI);
  while (phi < -M_PI)
    phi += (2 * M_PI);
  return phi;
}

// ------------ method called once each job just before starting event loop  ------------
void CastorClusterProducer::beginJob() { LogDebug("CastorClusterProducer") << "Starting CastorClusterProducer"; }

// ------------ method called once each job just after ending the event loop  ------------
void CastorClusterProducer::endJob() { LogDebug("CastorClusterProducer") << "Ending CastorClusterProducer"; }

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