JetProducers/src/CastorJetIDHelper.cc

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   // calculate Castor JetID properties
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     // loop over rechits
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     }  // end loop over rechits
0060
0061     nTowers_++;
0062   }
0063   //cout << "" << endl;
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 // help function to calculate phi within [-pi,+pi]
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 }