Castor/src/CastorClusterProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Castor
0004 // Class:      CastorClusterProducer
0005 //
0006 /**\class CastorClusterProducer CastorClusterProducer.cc RecoLocalCalo/Castor/src/CastorClusterProducer.cc
0007
0008  Description: CastorCluster Reconstruction Producer. Produces Clusters from Towers
0009  Implementation:
0010 */
0011
0012 //
0013 // Original Author:  Hans Van Haevermaet, Benoit Roland
0014 //         Created:  Wed Jul  9 14:00:40 CEST 2008
0015 //
0016 //
0017
0018 // system include
0019 #include <memory>
0020 #include <vector>
0021 #include <iostream>
0022 #include <TMath.h>
0023 #include <TRandom3.h>
0024
0025 // user include
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "DataFormats/Math/interface/Point3D.h"
0032
0033 // Castor Object include
0034 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
0035 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0036 #include "DataFormats/CastorReco/interface/CastorTower.h"
0037 #include "DataFormats/CastorReco/interface/CastorCluster.h"
0038 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0039 #include "DataFormats/Candidate/interface/Candidate.h"
0040 #include "DataFormats/JetReco/interface/BasicJet.h"
0041 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
0042 #include "DataFormats/JetReco/interface/Jet.h"
0043
0044 //
0045 // class decleration
0046 //
0047
0048 class CastorClusterProducer : public edm::global::EDProducer<> {
0049 public:
0050   explicit CastorClusterProducer(const edm::ParameterSet&);
0051   ~CastorClusterProducer() override;
0052
0053 private:
0054   void beginJob() override;
0055   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056   double phiangle(double testphi) const;
0057   void endJob() override;
0058
0059   // ----------member data ---------------------------
0060   typedef math::XYZPointD Point;
0061   typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
0062   typedef ROOT::Math::RhoZPhiPoint CellPoint;
0063   typedef std::vector<reco::CastorTower> CastorTowerCollection;
0064   typedef std::vector<reco::CastorCluster> CastorClusterCollection;
0065   std::string input_, basicjets_;
0066   edm::EDGetTokenT<CastorTowerCollection> tok_input_;
0067   edm::EDGetTokenT<reco::BasicJetCollection> tok_jets_;
0068   edm::EDGetTokenT<CastorTowerCollection> tok_tower_;
0069   bool clusteralgo_;
0070 };
0071
0072 //
0073 // constants, enums and typedefs
0074 //
0075
0076 //
0077 // static data member definitions
0078 //
0079
0080 //
0081 // constructor and destructor
0082 //
0083
0084 CastorClusterProducer::CastorClusterProducer(const edm::ParameterSet& iConfig)
0085     : input_(iConfig.getUntrackedParameter<std::string>("inputtowers", "")),
0086       basicjets_(iConfig.getUntrackedParameter<std::string>("basicjets", "")),
0087       clusteralgo_(iConfig.getUntrackedParameter<bool>("ClusterAlgo", false)) {
0088   // register for data access
0089   tok_input_ = consumes<CastorTowerCollection>(edm::InputTag(input_));
0090   tok_jets_ = consumes<reco::BasicJetCollection>(edm::InputTag(basicjets_));
0091   tok_tower_ = consumes<CastorTowerCollection>(edm::InputTag("CastorTowerReco"));
0092   // register your products
0093   produces<CastorClusterCollection>();
0094   // now do what ever other initialization is needed
0095 }
0096
0097 CastorClusterProducer::~CastorClusterProducer() {
0098   // do anything here that needs to be done at desctruction time
0099   // (e.g. close files, deallocate resources etc.)
0100 }
0101
0102 //
0103 // member functions
0104 //
0105
0106 // ------------ method called to produce the data  ------------
0107 void CastorClusterProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0108   using namespace edm;
0109   using namespace reco;
0110   using namespace TMath;
0111
0112   LogDebug("CastorClusterProducer") << "3. entering CastorClusterProducer";
0113
0114   if (!input_.empty()) {
0115     // Produce CastorClusters from CastorTowers
0116
0117     edm::Handle<CastorTowerCollection> InputTowers;
0118     iEvent.getByToken(tok_input_, InputTowers);
0119
0120     auto OutputClustersfromClusterAlgo = std::make_unique<CastorClusterCollection>();
0121
0122     // get and check input size
0123     int nTowers = InputTowers->size();
0124
0125     if (nTowers == 0)
0126       LogDebug("CastorClusterProducer") << "Warning: You are trying to run the Cluster algorithm with 0 input towers.";
0127
0128     CastorTowerRefVector posInputTowers, negInputTowers;
0129
0130     for (size_t i = 0; i < InputTowers->size(); ++i) {
0131       reco::CastorTowerRef tower_p = reco::CastorTowerRef(InputTowers, i);
0132       if (tower_p->eta() > 0.)
0133         posInputTowers.push_back(tower_p);
0134       if (tower_p->eta() < 0.)
0135         negInputTowers.push_back(tower_p);
0136     }
0137
0138     // build cluster from ClusterAlgo
0139     if (clusteralgo_ == true) {
0140       // code
0141       iEvent.put(std::move(OutputClustersfromClusterAlgo));
0142     }
0143   }
0144
0145   if (!basicjets_.empty()) {
0146     Handle<BasicJetCollection> bjCollection;
0147     iEvent.getByToken(tok_jets_, bjCollection);
0148
0149     Handle<CastorTowerCollection> ctCollection;
0150     iEvent.getByToken(tok_tower_, ctCollection);
0151
0152     auto OutputClustersfromBasicJets = std::make_unique<CastorClusterCollection>();
0153
0154     if (bjCollection->empty())
0155       LogDebug("CastorClusterProducer")
0156           << "Warning: You are trying to run the Cluster algorithm with 0 input basicjets.";
0157
0158     for (unsigned i = 0; i < bjCollection->size(); i++) {
0159       const BasicJet* bj = &(*bjCollection)[i];
0160
0161       double energy = bj->energy();
0162       TowerPoint temp(88.5, bj->eta(), bj->phi());
0163       Point position(temp);
0164       double emEnergy = 0.;
0165       double hadEnergy = 0.;
0166       double width = 0.;
0167       double depth = 0.;
0168       double fhot = 0.;
0169       double sigmaz = 0.;
0170       CastorTowerRefVector usedTowers;
0171       double zmean = 0.;
0172       double z2mean = 0.;
0173
0174       std::vector<CandidatePtr> ccp = bj->getJetConstituents();
0175       std::vector<CandidatePtr>::const_iterator itParticle;
0176       for (itParticle = ccp.begin(); itParticle != ccp.end(); ++itParticle) {
0177         const CastorTower* castorcand = dynamic_cast<const CastorTower*>(itParticle->get());
0178         //cout << " castortowercandidate reference energy = " << castorcand->castorTower()->energy() << endl;
0179         //cout << " castortowercandidate reference eta = " << castorcand->castorTower()->eta() << endl;
0180         //cout << " castortowercandidate reference phi = " << castorcand->castorTower()->phi() << endl;
0181         //cout << " castortowercandidate reference depth = " << castorcand->castorTower()->depth() << endl;
0182
0183         //CastorTowerCollection *ctc = new CastorTowerCollection();
0184         //ctc->push_back(*castorcand);
0185         //CastorTowerRef towerref = CastorTowerRef(ctc,0);
0186
0187         size_t thisone = 0;
0188         for (size_t l = 0; l < ctCollection->size(); l++) {
0189           const CastorTower ct = (*ctCollection)[l];
0190           if (std::abs(ct.phi() - castorcand->phi()) < 0.0001) {
0191             thisone = l;
0192           }
0193         }
0194
0195         CastorTowerRef towerref(ctCollection, thisone);
0196         usedTowers.push_back(towerref);
0197         emEnergy += castorcand->emEnergy();
0198         hadEnergy += castorcand->hadEnergy();
0199         depth += castorcand->depth() * castorcand->energy();
0200         width += pow(phiangle(castorcand->phi() - bj->phi()), 2) * castorcand->energy();
0201         fhot += castorcand->fhot() * castorcand->energy();
0202
0203         // loop over rechits
0204         for (edm::RefVector<edm::SortedCollection<CastorRecHit> >::iterator it = castorcand->rechitsBegin();
0205              it != castorcand->rechitsEnd();
0206              it++) {
0207           edm::Ref<edm::SortedCollection<CastorRecHit> > rechit_p = *it;
0208           double Erechit = rechit_p->energy();
0209           HcalCastorDetId id = rechit_p->id();
0210           int module = id.module();
0211           double zrechit = 0;
0212           if (module < 3)
0213             zrechit = -14390 - 24.75 - 49.5 * (module - 1);
0214           if (module > 2)
0215             zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
0216           zmean += Erechit * zrechit;
0217           z2mean += Erechit * zrechit * zrechit;
0218         }  // end loop over rechits
0219       }
0220       //cout << "" << endl;
0221
0222       depth = depth / energy;
0223       width = sqrt(width / energy);
0224       fhot = fhot / energy;
0225
0226       zmean = zmean / energy;
0227       z2mean = z2mean / energy;
0228       double sigmaz2 = z2mean - zmean * zmean;
0229       if (sigmaz2 > 0)
0230         sigmaz = sqrt(sigmaz2);
0231
0232       CastorCluster cc(
0233           energy, position, emEnergy, hadEnergy, emEnergy / energy, width, depth, fhot, sigmaz, usedTowers);
0234       OutputClustersfromBasicJets->push_back(cc);
0235     }
0236
0237     iEvent.put(std::move(OutputClustersfromBasicJets));
0238   }
0239 }
0240
0241 // help function to calculate phi within [-pi,+pi]
0242 double CastorClusterProducer::phiangle(double testphi) const {
0243   double phi = testphi;
0244   while (phi > M_PI)
0245     phi -= (2 * M_PI);
0246   while (phi < -M_PI)
0247     phi += (2 * M_PI);
0248   return phi;
0249 }
0250
0251 // ------------ method called once each job just before starting event loop  ------------
0252 void CastorClusterProducer::beginJob() { LogDebug("CastorClusterProducer") << "Starting CastorClusterProducer"; }
0253
0254 // ------------ method called once each job just after ending the event loop  ------------
0255 void CastorClusterProducer::endJob() { LogDebug("CastorClusterProducer") << "Ending CastorClusterProducer"; }
0256
0257 //define this as a plug-in
0258 DEFINE_FWK_MODULE(CastorClusterProducer);