/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoJets/JetProducers/src/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/CastorReco/interface/CastorCell.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 cells
                        for (CastorCell_iterator it = castorcand->cellsBegin(); it != castorcand->cellsEnd(); it++) {
                                CastorCellRef cell_p = *it;
                                math::XYZPointD rcell = cell_p->position();
                                double Ecell = cell_p->energy();
                                zmean += Ecell*cell_p->z();
                                z2mean += Ecell*cell_p->z()*cell_p->z();
                        } // end loop over cells
                        
                        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;
}