AlphaTVarProducer.cc

Go to the documentation of this file.

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Common/interface/Handle.h"

#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"

#include "DataFormats/METReco/interface/CaloMET.h"
#include "DataFormats/METReco/interface/CaloMETCollection.h"

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/Utilities/interface/InputTag.h"

#include "DQM/DataScouting/interface/AlphaTVarProducer.h"

#include "TVector3.h"

#include <memory>
#include <vector>

//
// constructors and destructor
//
AlphaTVarProducer::AlphaTVarProducer(const edm::ParameterSet& iConfig) :
  inputJetTag_    (iConfig.getParameter<edm::InputTag>("inputJetTag")){

  produces<std::vector<double> >();

  //set Token(-s)
  inputJetTagToken_ = consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("inputJetTag"));

  LogDebug("") << "Inputs: "
           << inputJetTag_.encode() << " ";
}

AlphaTVarProducer::~AlphaTVarProducer()
{
}

// ------------ method called to produce the data  ------------ 
void
AlphaTVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace std;
   using namespace edm;
   using namespace reco;


   // get hold of collection of objects
   edm::Handle<reco::CaloJetCollection> calojet_handle; 
   iEvent.getByToken(inputJetTagToken_, calojet_handle);

   std::auto_ptr<std::vector<double> > result(new std::vector<double>); 
   // check the the input collections are available
   if (calojet_handle.isValid()){
     std::vector<TLorentzVector> myJets;
     reco::CaloJetCollection::const_iterator jetIt;
     for(jetIt = calojet_handle->begin(); jetIt != calojet_handle->end(); jetIt++){
       TLorentzVector j; j.SetPtEtaPhiE(jetIt->pt(),jetIt->eta(), jetIt->phi(), jetIt->energy());
       myJets.push_back(j);
     }

     double alphaT = CalcAlphaT(myJets);
     double HT     = CalcHT(myJets);

     result->push_back(alphaT);
     result->push_back(HT);
   }
   iEvent.put(result);
}

double 
AlphaTVarProducer::CalcAlphaT(const std::vector<TLorentzVector>& jets){
    std::vector<double> ETs;
    TVector3 MHT{CalcMHT(jets),0.0,0.0};
    float HT = CalcHT(jets);
    //float HT = 0;
    for(unsigned int i = 0; i < jets.size(); i++){
        if(jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5) ETs.push_back(jets[i].Et());
        //HT += jets[i].Et();
    }
    if(ETs.size() < 2.) return 0.0;
    if(ETs.size() > 16.) return 0.0;
    float DHT = deltaHt(ETs);
    
    float AlphaT = alphaT(HT,DHT,MHT.Mag());
    
    return AlphaT;
    
}

double AlphaTVarProducer::deltaHt(const std::vector<double>& ETs) {
  if(ETs.size() > 16.) return 9999999;
    std::vector<double> diff( 1<<(ETs.size()-1) , 0. );
    for(unsigned i=0; i < diff.size(); i++)
        for(unsigned j=0; j < ETs.size(); j++)
            diff[i] += ETs[j] * ( 1 - 2 * (int(i>>j)&1) ) ;
    std::vector<double>::const_iterator it;
    double min=9999999;
    for(it = diff.begin(); it !=diff.end(); it++) if(*it<min) min = *it;
    return min;
}

double AlphaTVarProducer::alphaT(const double HT, const double DHT, const double MHT) {
    return 0.5 * ( HT - DHT ) / sqrt( HT*HT - MHT*MHT );
}

double 
AlphaTVarProducer::CalcHT(const std::vector<TLorentzVector>& jets){

    double HT = 0;
    for(unsigned int i = 0; i < jets.size(); i++){
        if(jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5) HT += jets[i].Et();
    }
    
    return HT;
    
}

double AlphaTVarProducer::CalcMHT(const std::vector<TLorentzVector>& jets){
    TVector3 MHT;
    for(unsigned int i = 0; i < jets.size(); i++){
        if(jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5) MHT -= jets[i].Vect();
    }
    MHT.SetZ(0.0);
    return MHT.Mag();
}

DEFINE_FWK_MODULE(AlphaTVarProducer);