PatExamples/plugins/MuonAnalyzer.cc

0001 /** \class ExampleMuonAnalyzer
0002  *  Analyzer of the muon objects
0003  *
0004  *  \author R. Bellan - CERN <riccardo.bellan@cern.ch>
0005  */
0006
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 #include "DataFormats/Math/interface/LorentzVector.h"
0009 #include "DataFormats/PatCandidates/interface/Muon.h"
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "TH1I.h"
0018 #include "TH1F.h"
0019 #include "TH2F.h"
0020
0021 class ExampleMuonAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0022 public:
0023   /// Constructor
0024   ExampleMuonAnalyzer(const edm::ParameterSet &pset);
0025
0026   /// Destructor
0027   ~ExampleMuonAnalyzer() override;
0028
0029   // Operations
0030
0031   void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override;
0032
0033   void beginJob() override;
0034   void endJob() override;
0035
0036 protected:
0037 private:
0038   edm::EDGetTokenT<pat::MuonCollection> theMuonToken;
0039
0040   // Histograms
0041   TH1I *hNMuons;
0042   TH1F *hPtRec;
0043   TH2F *hPtReso;
0044   TH1F *hEHcal;
0045
0046   TH1I *hMuonType;
0047   TH1F *hPtSTATKDiff;
0048
0049   // ID
0050   TH1F *hMuCaloCompatibility;
0051   TH1F *hMuSegCompatibility;
0052   TH1I *hChamberMatched;
0053   TH1I *hMuIdAlgo;
0054
0055   // Isolation
0056   TH1F *hMuIso03SumPt;
0057   TH1F *hMuIso03CaloComb;
0058
0059   TH1F *h4MuInvMass;
0060 };
0061
0062 using namespace std;
0063 using namespace edm;
0064
0065 /// Constructor
0066 ExampleMuonAnalyzer::ExampleMuonAnalyzer(const ParameterSet &pset) {
0067   theMuonToken = consumes<pat::MuonCollection>(pset.getParameter<InputTag>("MuonCollection"));
0068   usesResource(TFileService::kSharedResource);
0069 }
0070
0071 /// Destructor
0072 ExampleMuonAnalyzer::~ExampleMuonAnalyzer() {}
0073
0074 void ExampleMuonAnalyzer::beginJob() {
0075   // Book histograms
0076   edm::Service<TFileService> fileService;
0077   hPtRec = fileService->make<TH1F>("pT", "p_{T}^{rec}", 250, 0, 120);
0078   hPtReso = fileService->make<TH2F>("pT_Reso", "(p_{T}^{rec}-p_{T}^{sim})/p_{T}^{sim}", 250, 0, 120, 100, -0.2, 0.2);
0079   hNMuons = fileService->make<TH1I>("NMuons", "Number of muons per event", 20, 0, 20);
0080
0081   hEHcal = fileService->make<TH1F>("EHCal", "Energy deposit in HCAL", 100, 0, 10);
0082
0083   // ID
0084   hMuonType = fileService->make<TH1I>("MuonType", "Type of Muons", 5, 1, 6);
0085   hPtSTATKDiff = fileService->make<TH1F>("DiffPt_STA_TK", "p^{TK}_{T}-p^{STA}_T", 200, -50, 50);
0086   hMuCaloCompatibility = fileService->make<TH1F>("CaloCompatibility", "Muon HP using Calorimeters only", 100, 0, 1);
0087   hMuSegCompatibility = fileService->make<TH1F>("SegmentCompatibility", "Muon HP using segments only", 100, 0, 1);
0088   hChamberMatched = fileService->make<TH1I>("NumMatchedChamber", "Number of matched chambers", 7, 0, 7);
0089   hMuIdAlgo = fileService->make<TH1I>("MuonIDSelectors", "Results of muon id selectors", 13, 0, 13);
0090
0091   // Isolation
0092   hMuIso03SumPt = fileService->make<TH1F>("MuIso03SumPt", "Isolation #Delta(R)=0.3: SumPt", 200, 0, 10);
0093   hMuIso03CaloComb =
0094       fileService->make<TH1F>("MuIso03CaloComb", "Isolation #Delta(R)=0.3: 1.2*ECAL+0.8HCAL", 200, 0, 10);
0095
0096   // 4Mu invariant mass
0097   h4MuInvMass = fileService->make<TH1F>("InvMass4MuSystem", "Invariant mass of the 4 muons system", 200, 0, 500);
0098 }
0099
0100 void ExampleMuonAnalyzer::endJob() {}
0101
0102 void ExampleMuonAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0103   // Get the Muon collection
0104   Handle<pat::MuonCollection> muons;
0105   event.getByToken(theMuonToken, muons);
0106
0107   // How many muons in the event?
0108   hNMuons->Fill(muons->size());
0109
0110   pat::MuonCollection selectedMuons;
0111
0112   // Let's look inside the muon collection.
0113   for (pat::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0114     // pT spectra of muons
0115     hPtRec->Fill(muon->pt());
0116
0117     // what is the resolution in pt? Easy! We have the association with generated information
0118     //    cout<<muon->pt()<<" "<<muon->genParticle()->pt()<<endl;
0119     if (muon->genLepton() != nullptr) {
0120       double reso = (muon->pt() - muon->genLepton()->pt()) / muon->genLepton()->pt();
0121       hPtReso->Fill(muon->genLepton()->pt(), reso);
0122     }
0123
0124     // What is the energy deposit in HCal?
0125     if (muon->isEnergyValid())
0126       hEHcal->Fill(muon->calEnergy().had);
0127
0128     // Which type of muons in the collection?
0129     if (muon->isStandAloneMuon())
0130       if (muon->isGlobalMuon())
0131         if (muon->isTrackerMuon())
0132           hMuonType->Fill(1);  // STA + GLB + TM
0133         else
0134           hMuonType->Fill(2);  // STA + GLB
0135       else if (muon->isTrackerMuon())
0136         hMuonType->Fill(3);  // STA + TM
0137       else
0138         hMuonType->Fill(5);  // STA
0139     else if (muon->isTrackerMuon())
0140       hMuonType->Fill(4);  // TM
0141
0142     // ...mmm I want to study the relative resolution of the STA track with respect to the Tracker track.
0143     // or I want to look at a track stab
0144     if (muon->isGlobalMuon()) {
0145       double diff = muon->innerTrack()->pt() - muon->standAloneMuon()->pt();
0146       hPtSTATKDiff->Fill(diff);
0147     }
0148
0149     // Muon ID quantities
0150
0151     // Muon in CMS are usually MIP. What is the compatibility of a muon HP using only claorimeters?
0152     if (muon->isCaloCompatibilityValid())
0153       hMuCaloCompatibility->Fill(muon->caloCompatibility());
0154
0155     // The muon system can also be used just as only for ID. What is the compatibility of a muon HP using only calorimeters?
0156     hMuSegCompatibility->Fill(muon::segmentCompatibility(*muon));
0157
0158     // How many chambers have been associated to a muon track?
0159     hChamberMatched->Fill(muon->numberOfChambers());
0160     // If you look at MuonSegmentMatcher class you will see a lot of interesting quantities to look at!
0161     // you can get the list of matched info using matches()
0162
0163     // Muon ID selection. As described in AN-2008/098
0164     if (muon::isGoodMuon(*muon, muon::All))  // dummy options - always true
0165       hMuIdAlgo->Fill(0);
0166     if (muon::isGoodMuon(*muon, muon::AllStandAloneMuons))  // checks isStandAloneMuon flag
0167       hMuIdAlgo->Fill(1);
0168     if (muon::isGoodMuon(*muon, muon::AllTrackerMuons))  // checks isTrackerMuon flag
0169       hMuIdAlgo->Fill(2);
0170     if (muon::isGoodMuon(*muon, muon::TrackerMuonArbitrated))  // resolve ambiguity of sharing segments
0171       hMuIdAlgo->Fill(3);
0172     if (muon::isGoodMuon(*muon, muon::AllArbitrated))  // all muons with the tracker muon arbitrated
0173       hMuIdAlgo->Fill(4);
0174     if (muon::isGoodMuon(*muon, muon::GlobalMuonPromptTight))  // global muons with tighter fit requirements
0175       hMuIdAlgo->Fill(5);
0176     if (muon::isGoodMuon(*muon, muon::TMLastStationLoose))  // penetration depth loose selector
0177       hMuIdAlgo->Fill(6);
0178     if (muon::isGoodMuon(*muon, muon::TMLastStationTight))  // penetration depth tight selector
0179       hMuIdAlgo->Fill(7);
0180     if (muon::isGoodMuon(*muon, muon::TM2DCompatibilityLoose))  // likelihood based loose selector
0181       hMuIdAlgo->Fill(8);
0182     if (muon::isGoodMuon(*muon, muon::TM2DCompatibilityTight))  // likelihood based tight selector
0183       hMuIdAlgo->Fill(9);
0184     if (muon::isGoodMuon(*muon, muon::TMOneStationLoose))  // require one well matched segment
0185       hMuIdAlgo->Fill(10);
0186     if (muon::isGoodMuon(*muon, muon::TMOneStationTight))  // require one well matched segment
0187       hMuIdAlgo->Fill(11);
0188     if (muon::isGoodMuon(*muon,
0189                          muon::TMLastStationOptimizedLowPtLoose))  // combination of TMLastStation and TMOneStation
0190       hMuIdAlgo->Fill(12);
0191     if (muon::isGoodMuon(*muon,
0192                          muon::TMLastStationOptimizedLowPtTight))  // combination of TMLastStation and TMOneStation
0193       hMuIdAlgo->Fill(13);
0194
0195     // Isolation variables. There are many type of isolation. You can even build your own combining the output of
0196     // muon->isolationR03(). E.g.: 1.2*muon->isolationR03().emEt + 0.8*muon->isolationR03().hadEt
0197     // *** WARNING *** it is just an EXAMPLE!
0198     if (muon->isIsolationValid()) {
0199       hMuIso03CaloComb->Fill(1.2 * muon->isolationR03().emEt + 0.8 * muon->isolationR03().hadEt);
0200       hMuIso03SumPt->Fill(muon->isolationR03().sumPt);
0201     }
0202
0203     // OK, let see if we understood everything.
0204     // Suppose we are searching for H->ZZ->4mu.
0205     // In mean the 4 muons have/are:
0206     // high pt (but 1 out of 4 can be at quite low pt)
0207     // isolated
0208     // so, we can make some requirements
0209     if (muon->isolationR03().sumPt < 0.2) {
0210       if (muon->isGlobalMuon() || muon::isGoodMuon(*muon, muon::TM2DCompatibilityLoose) ||
0211           muon::isGoodMuon(*muon, muon::TMLastStationLoose))
0212         selectedMuons.push_back(*muon);
0213     }
0214   }
0215
0216   /// simple selection... Do not want to write here my super-secret Higgs analysis ;-)
0217   if (selectedMuons.size() == 4) {
0218     reco::Candidate::LorentzVector p4CM;
0219     for (pat::MuonCollection::const_iterator muon = selectedMuons.begin(); muon != selectedMuons.end(); ++muon) {
0220       p4CM = p4CM + muon->p4();
0221     }
0222     h4MuInvMass->Fill(p4CM.mass());
0223   }
0224 }
0225 DEFINE_FWK_MODULE(ExampleMuonAnalyzer);