JetAnalyzer_HeavyIons_matching.h

Go to the documentation of this file.

#ifndef JetAnalyzer_HeavyIons_matching_H
#define JetAnalyzer_HeavyIons_matching_H


//
// Jet Tester class for heavy ion jets. for DQM jet analysis monitoring 
// For CMSSW_7_5_X, especially reading background subtracted jets 
// author: Raghav Kunnawalkam Elayavalli,
//         April 6th 2015 
//         Rutgers University, email: raghav.k.e at CERN dot CH 
//
// The logic for the matching is taken from Pawan Kumar Netrakanti's macro analysis level macro available here
// https://github.com/pawannetrakanti/UserCode/blob/master/JetRAA/jetmatch.C 
//


#include <memory>
#include <fstream>
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"

#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/HLTReco/interface/TriggerObject.h"
#include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
#include "CommonTools/TriggerUtils/interface/GenericTriggerEventFlag.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "DataFormats/JetReco/interface/JPTJet.h"
#include "DataFormats/JetReco/interface/JPTJetCollection.h"
#include "DataFormats/JetReco/interface/PFJet.h"

#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Candidate/interface/CandidateWithRef.h"

// include the basic jet for the PuPF jets. 
#include "DataFormats/JetReco/interface/BasicJet.h"
#include "DataFormats/JetReco/interface/BasicJetCollection.h"
// include the pf candidates 
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
// include the voronoi subtraction
#include "DataFormats/HeavyIonEvent/interface/VoronoiBackground.h"
#include "RecoHI/HiJetAlgos/interface/UEParameters.h"
// include the centrality variables
#include "DataFormats/HeavyIonEvent/interface/Centrality.h"

#include "JetMETCorrections/Objects/interface/JetCorrector.h"
#include "RecoJets/JetProducers/interface/JetIDHelper.h"
#include "DQMOffline/JetMET/interface/JetMETDQMDCSFilter.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/DetId/interface/DetId.h"

#include "DataFormats/Scalers/interface/DcsStatus.h" 
#include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include <map>
#include <string>
#include <vector>
#include <algorithm>
#include <set>
#include <utility>


class MonitorElement;

class JetAnalyzer_HeavyIons_matching : public DQMEDAnalyzer {
 public:

  explicit JetAnalyzer_HeavyIons_matching (const edm::ParameterSet&);
  virtual ~JetAnalyzer_HeavyIons_matching();
  
  void analyze(const edm::Event&, const edm::EventSetup&) override; 
  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
  
  // calojet1token = vscalo 
  // calojet2token = pucalo 

 private:
  
  edm::InputTag   mInputJet1Collection;
  edm::InputTag   mInputJet2Collection;
  
  std::string     mOutputFile;
  std::string     JetType1;
  std::string     JetType2;
  double          mRecoJetPtThreshold;
  double          mRecoDelRMatch; 
  double          mRecoJetEtaCut; 
  
  //Tokens
  edm::EDGetTokenT<reco::CaloJetCollection> caloJet1Token_;
  edm::EDGetTokenT<reco::CaloJetCollection> caloJet2Token_;
  edm::EDGetTokenT<reco::PFJetCollection> pfJetsToken_;
  edm::EDGetTokenT<reco::BasicJetCollection> basicJetsToken_;
  edm::EDGetTokenT<reco::JPTJetCollection> jptJetsToken_;

  MonitorElement * mpT_ratio_Jet1Jet2;
  MonitorElement * mpT_Jet1_matched;
  MonitorElement * mpT_Jet2_matched;
  MonitorElement * mpT_Jet1_unmatched;
  MonitorElement * mpT_Jet2_unmatched; 

  MonitorElement * mHadEnergy_Jet1_unmatched;
  MonitorElement * mEmEnergy_Jet1_unmatched;
  MonitorElement * mChargedHadronEnergy_Jet1_unmatched;
  MonitorElement * mNeutralHadronEnergy_Jet1_unmatched;
  MonitorElement * mChargedEmEnergy_Jet1_unmatched;
  MonitorElement * mNeutralEmEnergy_Jet1_unmatched;
  MonitorElement * mChargedMuEnergy_Jet1_unmatched;

  MonitorElement * mHadEnergy_Jet2_unmatched;
  MonitorElement * mEmEnergy_Jet2_unmatched;
  MonitorElement * mChargedHadronEnergy_Jet2_unmatched;
  MonitorElement * mNeutralHadronEnergy_Jet2_unmatched;
  MonitorElement * mChargedEmEnergy_Jet2_unmatched;
  MonitorElement * mNeutralEmEnergy_Jet2_unmatched;
  MonitorElement * mChargedMuEnergy_Jet2_unmatched;

  struct MyJet{
    int   id;
    float pt;
    float eta;
    float phi;
  };


  typedef std::pair< MyJet, MyJet > ABJetPair;

  struct CompareMatchedJets {
    bool operator()(const ABJetPair &A1, const ABJetPair &A2){
      MyJet jet1_pair1 = A1.first; 
      MyJet jet2_pair1 = A1.second; 
      MyJet jet1_pair2 = A2.first; 
      MyJet jet2_pair2 = A2.second; 
      float delr1 = deltaRR(jet1_pair1.eta, jet1_pair1.phi, jet2_pair1.eta, jet2_pair1.phi);
      float delr2 = deltaRR(jet1_pair2.eta, jet1_pair2.phi, jet2_pair2.eta, jet2_pair2.phi);

      return ((delr1 < delr2) && (jet1_pair1.pt > jet1_pair2.pt));

    }
  };


  typedef std::multiset< ABJetPair, CompareMatchedJets > ABMatchedJets;
  typedef std::multiset< ABJetPair >::iterator ABItr;

  static float deltaRR(float eta1, float phi1, float eta2, float phi2)
  {
    float deta = eta1 - eta2;
    float dphi = fabs(phi1 - phi2);
    if(dphi > M_PI )dphi -= 2*M_PI;
    float dr = sqrt(deta*deta + dphi*dphi);
    return dr;
  }
  
  
  
};



#endif