CastorJetIDHelper.cc

Go to the documentation of this file.

#include "RecoJets/JetProducers/interface/CastorJetIDHelper.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/CastorReco/interface/CastorTower.h"
#include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
#include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"

#include "TMath.h"
#include <vector>
#include <numeric>
#include <iostream>


reco::helper::CastorJetIDHelper::CastorJetIDHelper()
{
 
  initValues();
}
  
void reco::helper::CastorJetIDHelper::initValues()
{
      emEnergy_ = 0.0;
      hadEnergy_ = 0.0; 
      fem_ = 0.0;
      width_ = 0.0;
      depth_ = 0.0;
      fhot_ = 0.0;
      sigmaz_ = 0.0;
      nTowers_ = 0;
}


void reco::helper::CastorJetIDHelper::calculate( const edm::Event& event, const reco::BasicJet &jet )
{
  initValues();
  
  // calculate Castor JetID properties
 
        double zmean = 0.;
        double z2mean = 0.;
    
        std::vector<CandidatePtr> ccp = jet.getJetConstituents();
        std::vector<CandidatePtr>::const_iterator itParticle;
        for (itParticle=ccp.begin();itParticle!=ccp.end();++itParticle){        
                const CastorTower* castorcand = dynamic_cast<const CastorTower*>(itParticle->get());
            emEnergy_ += castorcand->emEnergy();
            hadEnergy_ += castorcand->hadEnergy();
            depth_ += castorcand->depth()*castorcand->energy();
            width_ += pow(phiangle(castorcand->phi() - jet.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
            
            nTowers_++;
        }
        //cout << "" << endl;
        
        depth_ = depth_/jet.energy();
        width_ = sqrt(width_/jet.energy());
        fhot_ = fhot_/jet.energy();
        fem_ = emEnergy_/jet.energy();
        
        zmean = zmean/jet.energy();
            z2mean = z2mean/jet.energy();
            double sigmaz2 = z2mean - zmean*zmean;
            if(sigmaz2 > 0) sigmaz_ = sqrt(sigmaz2);

  
}

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