MCResolutions/src/MCResolutions.cc

// -*- C++ -*-
//
// Package:    MCResolutions
// Class:      MCResolutions
//
/**\class MCResolutions MCResolutions.cc JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
 
 Description: [one line class summary]
 
 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  Christian Sander,,,
//         Created:  Wed Oct  6 18:22:21 CEST 2010
// $Id: MCResolutions.cc,v 1.2 2010/10/07 11:38:49 csander Exp $
//
//

#include "JetResolutionFromMC/MCResolutions/interface/MCResolutions.h"

//
// constructors and destructor
//
MCResolutions::MCResolutions(const edm::ParameterSet& iConfig) {
    //now do what ever initialization is needed
    _jetTag = iConfig.getParameter<edm::InputTag> ("jetTag");
    _muonTag = iConfig.getParameter<edm::InputTag> ("muonTag");
    _genJetTag = iConfig.getParameter<edm::InputTag> ("genJetTag");
    _weightName = iConfig.getParameter<edm::InputTag> ("weightName");
    _bTag =  iConfig.getParameter<edm::InputTag> ("bTag");
    _EBRecHits = iConfig.getParameter<edm::InputTag> ("EBRecHits");
    _EERecHits = iConfig.getParameter<edm::InputTag> ("EERecHits");
    _deltaPhiDiJet = iConfig.getParameter<double> ("deltaPhiDiJet");
    _absCut3rdJet = iConfig.getParameter<double> ("absCut3rdJet");
    _relCut3rdJet = iConfig.getParameter<double> ("relCut3rdJet");
    _deltaRMatch = iConfig.getParameter<double> ("deltaRMatch");
    _deltaRMatchVeto = iConfig.getParameter<double> ("deltaRMatchVeto");
    _absPtVeto = iConfig.getParameter<double> ("absPtVeto");
    _relPtVeto = iConfig.getParameter<double> ("relPtVeto");
    _deltaRDeadECal = iConfig.getParameter<double> ("deltaRDeadECal");
    _GenJetPtCut = iConfig.getParameter<double> ("GenJetPtCut");
    _bTagCut = iConfig.getParameter<double> ("bTagCut");
    _bTagDeltaR = iConfig.getParameter<double> ("bTagDelatR");
    _maskedEcalChannelStatusThreshold = iConfig.getParameter<int> ("maskedEcalChannelStatusThreshold");

    EBinEdges.push_back(0);
    EBinEdges.push_back(20);
    EBinEdges.push_back(30);
    EBinEdges.push_back(50);
    EBinEdges.push_back(80);
    EBinEdges.push_back(120);
    EBinEdges.push_back(170);
    EBinEdges.push_back(230);
    EBinEdges.push_back(300);
    EBinEdges.push_back(380);
    EBinEdges.push_back(470);
    EBinEdges.push_back(570);
    EBinEdges.push_back(680);
    EBinEdges.push_back(800);
    EBinEdges.push_back(1000);
    EBinEdges.push_back(1300);
    EBinEdges.push_back(1700);
    EBinEdges.push_back(2200);

    PtBinEdges.push_back(0);
    PtBinEdges.push_back(20);
    PtBinEdges.push_back(30);
    PtBinEdges.push_back(50);
    PtBinEdges.push_back(80);
    PtBinEdges.push_back(120);
    PtBinEdges.push_back(170);
    PtBinEdges.push_back(230);
    PtBinEdges.push_back(300);
    PtBinEdges.push_back(380);
    PtBinEdges.push_back(470);
    PtBinEdges.push_back(570);
    PtBinEdges.push_back(680);
    PtBinEdges.push_back(800);
    PtBinEdges.push_back(1000);
    PtBinEdges.push_back(1300);
    PtBinEdges.push_back(1700);
    PtBinEdges.push_back(2200);

    EtaBinEdges.push_back(0.0);
    EtaBinEdges.push_back(1.4);
    EtaBinEdges.push_back(2.6);
    EtaBinEdges.push_back(3.2);
    EtaBinEdges.push_back(4.5);
    EtaBinEdges.push_back(5.0);

    h_tot_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
    h_b_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
    h_dead_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
    for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_JetResPt_Pt.begin(); it != h_tot_JetResPt_Pt.end(); ++it) {
        it->resize(PtBinEdges.size() - 1);
    }
    for (std::vector<std::vector<TH1F*> >::iterator it = h_b_JetResPt_Pt.begin(); it != h_b_JetResPt_Pt.end(); ++it) {
        it->resize(PtBinEdges.size() - 1);
    }
    for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_JetResPt_Pt.begin(); it != h_dead_JetResPt_Pt.end(); ++it) {
        it->resize(PtBinEdges.size() - 1);
    }


    h_tot_JetResPt_E.resize(EtaBinEdges.size() - 1);
    h_b_JetResPt_E.resize(EtaBinEdges.size() - 1);
    h_dead_JetResPt_E.resize(EtaBinEdges.size() - 1);
    for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_JetResPt_E.begin(); it != h_tot_JetResPt_E.end(); ++it) {
        it->resize(EBinEdges.size() - 1);
    }
    for (std::vector<std::vector<TH1F*> >::iterator it = h_b_JetResPt_E.begin(); it != h_b_JetResPt_E.end(); ++it) {
        it->resize(EBinEdges.size() - 1);
    }
    for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_JetResPt_E.begin(); it != h_dead_JetResPt_E.end(); ++it) {
        it->resize(EBinEdges.size() - 1);
    }

}


MCResolutions::~MCResolutions() {

    // do anything here that needs to be done at desctruction time
    // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
MCResolutions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
    using namespace std;

    // Event setup
    envSet(iSetup);

    getChannelStatusMaps();

    //Weight
    edm::Handle<double> event_weight;
    iEvent.getByLabel(_weightName, event_weight);
    weight = (event_weight.isValid() ? (*event_weight) : 1.0);

    //GenJets
    edm::Handle<edm::View<reco::GenJet> > Jets_gen;
    iEvent.getByLabel(_genJetTag, Jets_gen);

    //PATJets
    edm::Handle<edm::View<reco::CaloJet> > Jets_rec;
    iEvent.getByLabel(_jetTag, Jets_rec);

    //softMuonBJetTags
    edm::Handle<reco::JetTagCollection> bTagHandle;
    iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle);
    const reco::JetTagCollection & bTags = *(bTagHandle.product());

    for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->
            begin();
            it != Jets_gen->end();
            ++it) {
        if (it->pt()<_GenJetPtCut)
            continue;
        const reco::CaloJet* matchedJet = 0;
        const reco::CaloJet* nearestJet = 0;
        double dRmatched = 999.;
        double dRnearest = 999.;
        for (edm::View<reco::CaloJet>::const_iterator jt = Jets_rec->
                begin();
                jt != Jets_rec->end();
                ++jt) {
            double dR = deltaR(*it, *jt);
            if (dR < dRmatched) {
                nearestJet = matchedJet;
                dRnearest = dRmatched;
                matchedJet = &(*jt);
                dRmatched = dR;
            } else if (dR < dRnearest) {
                nearestJet = &(*jt);
                dRnearest = dR;
            }
        }
        //// look if there is no further significant CaloJet near the genJet
        if (dRmatched < _deltaRMatch && ( nearestJet == 0 || dRnearest > _deltaRMatchVeto || (nearestJet->pt()
                                          < _absPtVeto && nearestJet->pt()/matchedJet->pt() < _relPtVeto)
                                        )
           ) {
            //cout << dRmatched << endl;
            double res = matchedJet->pt() / it->pt();
            h_tot_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
            h_tot_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
            //// look if matched jet can be matched to a b-tagged jet
            for (int i = 0; i != (int) bTags.size(); ++i) {
                double dR_btag = deltaR(*matchedJet, *(bTags[i].first));
                if (dR_btag < _bTagDeltaR && bTags[i].second > _bTagCut) {
                    //cout << "matched jet has b tag discriminator = "<<bTags[i].second << endl;
                    h_b_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
                    h_b_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
                    break;
                }
            }
            //// look if matched jet points into direction of masked ECAL cluster
            for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt != EcalAllDeadChannelsValMap.end(); ++kt) {
                math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
                double dR_dead = deltaR(*matchedJet, Evec);
                if (dR_dead < _deltaRDeadECal) {
                    //cout << "jet is pointing to dead Ecal cluster (eta, phi)= "<< it->eta() << ", " << it->phi() << endl;
                    h_dead_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
                    h_dead_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
                    break;
                }
            }
        }
    }

}


// ------------ method called once each job just before starting event loop  ------------
void
MCResolutions::beginJob() {
    edm::Service<TFileService> fs;
    if (!fs) {
        throw edm::Exception(edm::errors::Configuration, "TFile Service is not registered in cfg file");
    }

    for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
        for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
            char hname[100];
            sprintf(hname, "h_tot_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
            h_tot_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
            h_tot_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
            sprintf(hname, "h_b_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
            h_b_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
            h_b_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
            sprintf(hname, "h_dead_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
            h_dead_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
            h_dead_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
        }
    }

    for (unsigned int i_e = 0; i_e < EBinEdges.size() - 1; ++i_e) {
        for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
            char hname[100];
            sprintf(hname, "h_tot_ResponsePt_E%i_Eta%i", i_e, i_eta);
            h_tot_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
            h_tot_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
            sprintf(hname, "h_b_ResponsePt_E%i_Eta%i", i_e, i_eta);
            h_b_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
            h_b_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
            sprintf(hname, "h_dead_ResponsePt_E%i_Eta%i", i_e, i_eta);
            h_dead_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
            h_dead_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
        }
    }

    mapsReady = false;
}

// ------------ method called once each job just after ending the event loop  ------------
void
MCResolutions::endJob() {}

int MCResolutions::EBin(const double& e) {
    int i_e = -1;
    for (std::vector<double>::const_iterator it = EBinEdges.begin(); it != EBinEdges.end(); ++it) {
        if ((*it) > e)
            break;
        ++i_e;
    }
    if (i_e < 0)
        i_e = 0;
    if (i_e > (int)EBinEdges.size() - 2)
        i_e = (int)EBinEdges.size() - 2;

    return i_e;
}

int MCResolutions::PtBin(const double& pt) {
    int i_pt = -1;
    for (std::vector<double>::const_iterator it = PtBinEdges.begin(); it != PtBinEdges.end(); ++it) {
        if ((*it) > pt)
            break;
        ++i_pt;
    }
    if (i_pt < 0)
        i_pt = 0;
    if (i_pt > (int)PtBinEdges.size() - 2)
        i_pt = (int)PtBinEdges.size() - 2;

    return i_pt;
}

int MCResolutions::EtaBin(const double& eta) {
    int i_eta = -1;
    for (std::vector<double>::const_iterator it = EtaBinEdges.begin(); it != EtaBinEdges.end(); ++it) {
        if ((*it) > fabs(eta))
            break;
        ++i_eta;
    }
    if (i_eta < 0)
        i_eta = 0;
    if (i_eta > (int)EtaBinEdges.size() - 2)
        i_eta = (int)EtaBinEdges.size() - 2;
    return i_eta;
}

void MCResolutions::envSet(const edm::EventSetup& iSetup) {

    ecalScale.setEventSetup( iSetup );

    iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
    iSetup.get<CaloGeometryRecord> ().get(geometry);

    if ( !ecalStatus.isValid() )
        throw "Failed to get ECAL channel status!";
    if ( !geometry.isValid() )
        throw "Failed to get the geometry!";

}

int MCResolutions::getChannelStatusMaps() {

    if( mapsReady )
        return -1;

    EcalAllDeadChannelsValMap.clear();

    // Loop over EB ...
    for( int ieta=-85; ieta<=85; ieta++ ) {
        for( int iphi=0; iphi<=360; iphi++ ) {
            if(! EBDetId::validDetId( ieta, iphi ) )
                continue;

            const EBDetId detid = EBDetId( ieta, iphi, EBDetId::ETAPHIMODE );
            EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
            // refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
            int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;
            //std::cout << ieta << " " << iphi << " " << status << std:: endl;
            const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
            const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
            double eta = cellGeom->getPosition ().eta ();
            double phi = cellGeom->getPosition ().phi ();
            double theta = cellGeom->getPosition().theta();

            if(status >= _maskedEcalChannelStatusThreshold) {
                std::vector<double> valVec;
                valVec.push_back(eta);
                valVec.push_back(phi);
                valVec.push_back(theta);
                //std::cout << eta << " " << phi << std::endl;
                EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
            }
        } // end loop iphi
    } // end loop ieta

    // Loop over EE detid
    for( int ix=0; ix<=100; ix++ ) {
        for( int iy=0; iy<=100; iy++ ) {
            for( int iz=-1; iz<=1; iz++ ) {
                if(iz==0)
                    continue;
                if(! EEDetId::validDetId( ix, iy, iz ) )
                    continue;

                const EEDetId detid = EEDetId( ix, iy, iz, EEDetId::XYMODE );
                EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
                int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;

                const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
                const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
                double eta = cellGeom->getPosition ().eta () ;
                double phi = cellGeom->getPosition ().phi () ;
                double theta = cellGeom->getPosition().theta();

                if(status >= _maskedEcalChannelStatusThreshold) {
                    std::vector<double> valVec;
                    valVec.push_back(eta);
                    valVec.push_back(phi);
                    valVec.push_back(theta);
                    //std::cout << eta << " " << phi << std::endl;
                    EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
                }
            } // end loop iz
        } // end loop iy
    } // end loop ix

    mapsReady = true;
    return 1;
}

//define this as a plug-in
DEFINE_FWK_MODULE(MCResolutions);
Revision:	1.3
Committed:	Fri Oct 8 17:48:45 2010 UTC (14 years, 6 months ago) by csander
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+319 -15 lines
Log Message:	first running version
#	User	Rev	Content
1	csander	1.1	// -- C++ --
2			//
3			// Package: MCResolutions
4			// Class: MCResolutions
5	csander	1.2	//
6	csander	1.1	/**\class MCResolutions MCResolutions.cc JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
7	csander	1.2
8	csander	1.1	Description: [one line class summary]
9	csander	1.2
10	csander	1.1	Implementation:
11			[Notes on implementation]
12			*/
13			//
14			// Original Author: Christian Sander,,,
15			// Created: Wed Oct 6 18:22:21 CEST 2010
16	csander	1.3	// $Id: MCResolutions.cc,v 1.2 2010/10/07 11:38:49 csander Exp $
17	csander	1.1	//
18			//
19
20	csander	1.2	#include "JetResolutionFromMC/MCResolutions/interface/MCResolutions.h"
21	csander	1.1
22			//
23			// constructors and destructor
24			//
25	csander	1.3	MCResolutions::MCResolutions(const edm::ParameterSet& iConfig) {
26	csander	1.2	//now do what ever initialization is needed
27			_jetTag = iConfig.getParameter<edm::InputTag> ("jetTag");
28			_muonTag = iConfig.getParameter<edm::InputTag> ("muonTag");
29			_genJetTag = iConfig.getParameter<edm::InputTag> ("genJetTag");
30			_weightName = iConfig.getParameter<edm::InputTag> ("weightName");
31	csander	1.3	_bTag = iConfig.getParameter<edm::InputTag> ("bTag");
32	csander	1.2	_EBRecHits = iConfig.getParameter<edm::InputTag> ("EBRecHits");
33			_EERecHits = iConfig.getParameter<edm::InputTag> ("EERecHits");
34			_deltaPhiDiJet = iConfig.getParameter<double> ("deltaPhiDiJet");
35			_absCut3rdJet = iConfig.getParameter<double> ("absCut3rdJet");
36			_relCut3rdJet = iConfig.getParameter<double> ("relCut3rdJet");
37			_deltaRMatch = iConfig.getParameter<double> ("deltaRMatch");
38			_deltaRMatchVeto = iConfig.getParameter<double> ("deltaRMatchVeto");
39			_absPtVeto = iConfig.getParameter<double> ("absPtVeto");
40			_relPtVeto = iConfig.getParameter<double> ("relPtVeto");
41			_deltaRDeadECal = iConfig.getParameter<double> ("deltaRDeadECal");
42	csander	1.3	_GenJetPtCut = iConfig.getParameter<double> ("GenJetPtCut");
43			_bTagCut = iConfig.getParameter<double> ("bTagCut");
44			_bTagDeltaR = iConfig.getParameter<double> ("bTagDelatR");
45			_maskedEcalChannelStatusThreshold = iConfig.getParameter<int> ("maskedEcalChannelStatusThreshold");
46
47			EBinEdges.push_back(0);
48			EBinEdges.push_back(20);
49			EBinEdges.push_back(30);
50			EBinEdges.push_back(50);
51			EBinEdges.push_back(80);
52			EBinEdges.push_back(120);
53			EBinEdges.push_back(170);
54			EBinEdges.push_back(230);
55			EBinEdges.push_back(300);
56			EBinEdges.push_back(380);
57			EBinEdges.push_back(470);
58			EBinEdges.push_back(570);
59			EBinEdges.push_back(680);
60			EBinEdges.push_back(800);
61			EBinEdges.push_back(1000);
62			EBinEdges.push_back(1300);
63			EBinEdges.push_back(1700);
64			EBinEdges.push_back(2200);
65
66			PtBinEdges.push_back(0);
67			PtBinEdges.push_back(20);
68			PtBinEdges.push_back(30);
69			PtBinEdges.push_back(50);
70			PtBinEdges.push_back(80);
71			PtBinEdges.push_back(120);
72			PtBinEdges.push_back(170);
73			PtBinEdges.push_back(230);
74			PtBinEdges.push_back(300);
75			PtBinEdges.push_back(380);
76			PtBinEdges.push_back(470);
77			PtBinEdges.push_back(570);
78			PtBinEdges.push_back(680);
79			PtBinEdges.push_back(800);
80			PtBinEdges.push_back(1000);
81			PtBinEdges.push_back(1300);
82			PtBinEdges.push_back(1700);
83			PtBinEdges.push_back(2200);
84
85			EtaBinEdges.push_back(0.0);
86			EtaBinEdges.push_back(1.4);
87			EtaBinEdges.push_back(2.6);
88			EtaBinEdges.push_back(3.2);
89			EtaBinEdges.push_back(4.5);
90			EtaBinEdges.push_back(5.0);
91
92			h_tot_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
93			h_b_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
94			h_dead_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
95			for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_JetResPt_Pt.begin(); it != h_tot_JetResPt_Pt.end(); ++it) {
96			it->resize(PtBinEdges.size() - 1);
97			}
98			for (std::vector<std::vector<TH1F*> >::iterator it = h_b_JetResPt_Pt.begin(); it != h_b_JetResPt_Pt.end(); ++it) {
99			it->resize(PtBinEdges.size() - 1);
100			}
101			for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_JetResPt_Pt.begin(); it != h_dead_JetResPt_Pt.end(); ++it) {
102			it->resize(PtBinEdges.size() - 1);
103			}
104
105
106			h_tot_JetResPt_E.resize(EtaBinEdges.size() - 1);
107			h_b_JetResPt_E.resize(EtaBinEdges.size() - 1);
108			h_dead_JetResPt_E.resize(EtaBinEdges.size() - 1);
109			for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_JetResPt_E.begin(); it != h_tot_JetResPt_E.end(); ++it) {
110			it->resize(EBinEdges.size() - 1);
111			}
112			for (std::vector<std::vector<TH1F*> >::iterator it = h_b_JetResPt_E.begin(); it != h_b_JetResPt_E.end(); ++it) {
113			it->resize(EBinEdges.size() - 1);
114			}
115			for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_JetResPt_E.begin(); it != h_dead_JetResPt_E.end(); ++it) {
116			it->resize(EBinEdges.size() - 1);
117			}
118	csander	1.1
119			}
120
121
122	csander	1.3	MCResolutions::~MCResolutions() {
123	csander	1.2
124			// do anything here that needs to be done at desctruction time
125			// (e.g. close files, deallocate resources etc.)
126	csander	1.1
127			}
128
129
130			//
131			// member functions
132			//
133
134			// ------------ method called to for each event ------------
135			void
136	csander	1.3	MCResolutions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
137			using namespace std;
138
139			// Event setup
140			envSet(iSetup);
141
142			getChannelStatusMaps();
143	csander	1.1
144	csander	1.2	//Weight
145			edm::Handle<double> event_weight;
146			iEvent.getByLabel(_weightName, event_weight);
147			weight = (event_weight.isValid() ? (*event_weight) : 1.0);
148
149			//GenJets
150			edm::Handle<edm::View<reco::GenJet> > Jets_gen;
151	csander	1.3	iEvent.getByLabel(_genJetTag, Jets_gen);
152	csander	1.2
153			//PATJets
154			edm::Handle<edm::View<reco::CaloJet> > Jets_rec;
155	csander	1.3	iEvent.getByLabel(_jetTag, Jets_rec);
156
157			//softMuonBJetTags
158			edm::Handle<reco::JetTagCollection> bTagHandle;
159			iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle);
160			const reco::JetTagCollection & bTags = *(bTagHandle.product());
161
162			for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->
163			begin();
164			it != Jets_gen->end();
165			++it) {
166			if (it->pt()<_GenJetPtCut)
167			continue;
168			const reco::CaloJet* matchedJet = 0;
169			const reco::CaloJet* nearestJet = 0;
170			double dRmatched = 999.;
171			double dRnearest = 999.;
172			for (edm::View<reco::CaloJet>::const_iterator jt = Jets_rec->
173			begin();
174			jt != Jets_rec->end();
175			++jt) {
176			double dR = deltaR(it, jt);
177			if (dR < dRmatched) {
178			nearestJet = matchedJet;
179			dRnearest = dRmatched;
180			matchedJet = &(*jt);
181			dRmatched = dR;
182			} else if (dR < dRnearest) {
183			nearestJet = &(*jt);
184			dRnearest = dR;
185			}
186			}
187			//// look if there is no further significant CaloJet near the genJet
188			if (dRmatched < _deltaRMatch && ( nearestJet == 0 \|\| dRnearest > _deltaRMatchVeto \|\| (nearestJet->pt()
189			< _absPtVeto && nearestJet->pt()/matchedJet->pt() < _relPtVeto)
190			)
191			) {
192			//cout << dRmatched << endl;
193			double res = matchedJet->pt() / it->pt();
194			h_tot_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
195			h_tot_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
196			//// look if matched jet can be matched to a b-tagged jet
197			for (int i = 0; i != (int) bTags.size(); ++i) {
198			double dR_btag = deltaR(matchedJet, (bTags[i].first));
199			if (dR_btag < _bTagDeltaR && bTags[i].second > _bTagCut) {
200			//cout << "matched jet has b tag discriminator = "<<bTags[i].second << endl;
201			h_b_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
202			h_b_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
203			break;
204			}
205			}
206			//// look if matched jet points into direction of masked ECAL cluster
207			for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt != EcalAllDeadChannelsValMap.end(); ++kt) {
208			math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
209			double dR_dead = deltaR(*matchedJet, Evec);
210			if (dR_dead < _deltaRDeadECal) {
211			//cout << "jet is pointing to dead Ecal cluster (eta, phi)= "<< it->eta() << ", " << it->phi() << endl;
212			h_dead_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res);
213			h_dead_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res);
214			break;
215			}
216			}
217			}
218			}
219	csander	1.1
220			}
221
222
223			// ------------ method called once each job just before starting event loop ------------
224	csander	1.2	void
225	csander	1.3	MCResolutions::beginJob() {
226			edm::Service<TFileService> fs;
227			if (!fs) {
228			throw edm::Exception(edm::errors::Configuration, "TFile Service is not registered in cfg file");
229			}
230
231			for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
232			for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
233			char hname[100];
234			sprintf(hname, "h_tot_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
235			h_tot_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
236			h_tot_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
237			sprintf(hname, "h_b_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
238			h_b_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
239			h_b_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
240			sprintf(hname, "h_dead_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
241			h_dead_JetResPt_Pt.at(i_eta).at(i_pt) = fs->make<TH1F> (hname, hname, 300, 0., 3.);
242			h_dead_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
243			}
244			}
245
246			for (unsigned int i_e = 0; i_e < EBinEdges.size() - 1; ++i_e) {
247			for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
248			char hname[100];
249			sprintf(hname, "h_tot_ResponsePt_E%i_Eta%i", i_e, i_eta);
250			h_tot_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
251			h_tot_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
252			sprintf(hname, "h_b_ResponsePt_E%i_Eta%i", i_e, i_eta);
253			h_b_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
254			h_b_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
255			sprintf(hname, "h_dead_ResponsePt_E%i_Eta%i", i_e, i_eta);
256			h_dead_JetResPt_E.at(i_eta).at(i_e) = fs->make<TH1F>(hname, hname, 300, 0., 3.);
257			h_dead_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
258			}
259			}
260
261			mapsReady = false;
262			}
263	csander	1.1
264			// ------------ method called once each job just after ending the event loop ------------
265	csander	1.2	void
266	csander	1.3	MCResolutions::endJob() {}
267
268			int MCResolutions::EBin(const double& e) {
269			int i_e = -1;
270			for (std::vector<double>::const_iterator it = EBinEdges.begin(); it != EBinEdges.end(); ++it) {
271			if ((*it) > e)
272			break;
273			++i_e;
274			}
275			if (i_e < 0)
276			i_e = 0;
277			if (i_e > (int)EBinEdges.size() - 2)
278			i_e = (int)EBinEdges.size() - 2;
279
280			return i_e;
281			}
282
283			int MCResolutions::PtBin(const double& pt) {
284			int i_pt = -1;
285			for (std::vector<double>::const_iterator it = PtBinEdges.begin(); it != PtBinEdges.end(); ++it) {
286			if ((*it) > pt)
287			break;
288			++i_pt;
289			}
290			if (i_pt < 0)
291			i_pt = 0;
292			if (i_pt > (int)PtBinEdges.size() - 2)
293			i_pt = (int)PtBinEdges.size() - 2;
294
295			return i_pt;
296			}
297
298			int MCResolutions::EtaBin(const double& eta) {
299			int i_eta = -1;
300			for (std::vector<double>::const_iterator it = EtaBinEdges.begin(); it != EtaBinEdges.end(); ++it) {
301			if ((*it) > fabs(eta))
302			break;
303			++i_eta;
304			}
305			if (i_eta < 0)
306			i_eta = 0;
307			if (i_eta > (int)EtaBinEdges.size() - 2)
308			i_eta = (int)EtaBinEdges.size() - 2;
309			return i_eta;
310			}
311
312			void MCResolutions::envSet(const edm::EventSetup& iSetup) {
313
314			ecalScale.setEventSetup( iSetup );
315
316			iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
317			iSetup.get<CaloGeometryRecord> ().get(geometry);
318
319			if ( !ecalStatus.isValid() )
320			throw "Failed to get ECAL channel status!";
321			if ( !geometry.isValid() )
322			throw "Failed to get the geometry!";
323
324			}
325
326			int MCResolutions::getChannelStatusMaps() {
327
328			if( mapsReady )
329			return -1;
330
331			EcalAllDeadChannelsValMap.clear();
332
333			// Loop over EB ...
334			for( int ieta=-85; ieta<=85; ieta++ ) {
335			for( int iphi=0; iphi<=360; iphi++ ) {
336			if(! EBDetId::validDetId( ieta, iphi ) )
337			continue;
338
339			const EBDetId detid = EBDetId( ieta, iphi, EBDetId::ETAPHIMODE );
340			EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
341			// refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
342			int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;
343			//std::cout << ieta << " " << iphi << " " << status << std:: endl;
344			const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
345			const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
346			double eta = cellGeom->getPosition ().eta ();
347			double phi = cellGeom->getPosition ().phi ();
348			double theta = cellGeom->getPosition().theta();
349
350			if(status >= _maskedEcalChannelStatusThreshold) {
351			std::vector<double> valVec;
352			valVec.push_back(eta);
353			valVec.push_back(phi);
354			valVec.push_back(theta);
355			//std::cout << eta << " " << phi << std::endl;
356			EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
357			}
358			} // end loop iphi
359			} // end loop ieta
360
361			// Loop over EE detid
362			for( int ix=0; ix<=100; ix++ ) {
363			for( int iy=0; iy<=100; iy++ ) {
364			for( int iz=-1; iz<=1; iz++ ) {
365			if(iz==0)
366			continue;
367			if(! EEDetId::validDetId( ix, iy, iz ) )
368			continue;
369
370			const EEDetId detid = EEDetId( ix, iy, iz, EEDetId::XYMODE );
371			EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
372			int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;
373
374			const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
375			const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
376			double eta = cellGeom->getPosition ().eta () ;
377			double phi = cellGeom->getPosition ().phi () ;
378			double theta = cellGeom->getPosition().theta();
379
380			if(status >= _maskedEcalChannelStatusThreshold) {
381			std::vector<double> valVec;
382			valVec.push_back(eta);
383			valVec.push_back(phi);
384			valVec.push_back(theta);
385			//std::cout << eta << " " << phi << std::endl;
386			EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
387			}
388			} // end loop iz
389			} // end loop iy
390			} // end loop ix
391
392			mapsReady = true;
393			return 1;
394			}
395	csander	1.1
396			//define this as a plug-in
397			DEFINE_FWK_MODULE(MCResolutions);