/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQM/DataScouting/src/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> >();

  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.getByLabel(inputJetTag_,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(std::vector<TLorentzVector> jets){
        std::vector<double> ETs;
        TVector3 MHT = CalcMHT(jets);
        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(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(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);