File indexing completed on 2023-03-17 11:18:34
0001 #include "RecoJets/JetProducers/interface/CastorJetIDHelper.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/CastorReco/interface/CastorTower.h"
0005 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
0006 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0007
0008 #include "TMath.h"
0009 #include <vector>
0010 #include <numeric>
0011 #include <iostream>
0012
0013 reco::helper::CastorJetIDHelper::CastorJetIDHelper() { initValues(); }
0014
0015 void reco::helper::CastorJetIDHelper::initValues() {
0016 emEnergy_ = 0.0;
0017 hadEnergy_ = 0.0;
0018 fem_ = 0.0;
0019 width_ = 0.0;
0020 depth_ = 0.0;
0021 fhot_ = 0.0;
0022 sigmaz_ = 0.0;
0023 nTowers_ = 0;
0024 }
0025
0026 void reco::helper::CastorJetIDHelper::calculate(const edm::Event& event, const reco::BasicJet& jet) {
0027 initValues();
0028
0029
0030
0031 double zmean = 0.;
0032 double z2mean = 0.;
0033
0034 std::vector<CandidatePtr> ccp = jet.getJetConstituents();
0035 std::vector<CandidatePtr>::const_iterator itParticle;
0036 for (itParticle = ccp.begin(); itParticle != ccp.end(); ++itParticle) {
0037 const CastorTower* castorcand = dynamic_cast<const CastorTower*>(itParticle->get());
0038 emEnergy_ += castorcand->emEnergy();
0039 hadEnergy_ += castorcand->hadEnergy();
0040 depth_ += castorcand->depth() * castorcand->energy();
0041 width_ += pow(phiangle(castorcand->phi() - jet.phi()), 2) * castorcand->energy();
0042 fhot_ += castorcand->fhot() * castorcand->energy();
0043
0044
0045 for (edm::RefVector<edm::SortedCollection<CastorRecHit> >::iterator it = castorcand->rechitsBegin();
0046 it != castorcand->rechitsEnd();
0047 it++) {
0048 edm::Ref<edm::SortedCollection<CastorRecHit> > rechit_p = *it;
0049 double Erechit = rechit_p->energy();
0050 HcalCastorDetId id = rechit_p->id();
0051 int module = id.module();
0052 double zrechit = 0;
0053 if (module < 3)
0054 zrechit = -14390 - 24.75 - 49.5 * (module - 1);
0055 if (module > 2)
0056 zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
0057 zmean += Erechit * zrechit;
0058 z2mean += Erechit * zrechit * zrechit;
0059 }
0060
0061 nTowers_++;
0062 }
0063
0064
0065 depth_ = depth_ / jet.energy();
0066 width_ = sqrt(width_ / jet.energy());
0067 fhot_ = fhot_ / jet.energy();
0068 fem_ = emEnergy_ / jet.energy();
0069
0070 zmean = zmean / jet.energy();
0071 z2mean = z2mean / jet.energy();
0072 double sigmaz2 = z2mean - zmean * zmean;
0073 if (sigmaz2 > 0)
0074 sigmaz_ = sqrt(sigmaz2);
0075 }
0076
0077
0078 double reco::helper::CastorJetIDHelper::phiangle(double testphi) {
0079 double phi = testphi;
0080 while (phi > M_PI)
0081 phi -= (2 * M_PI);
0082 while (phi < -M_PI)
0083 phi += (2 * M_PI);
0084 return phi;
0085 }