CMSSW Reference Manual

00001 /* \class ZMuMuEfficiency
00002  * 
00003  * author: Pasquale Noli
00004  * revised by Salvatore di Guida
00005  *
00006  * Efficiency of reconstruction tracker and muon Chamber
00007  *
00008  */
00009 
00010 #include "DataFormats/Common/interface/AssociationVector.h"
00011 #include "FWCore/Framework/interface/EDAnalyzer.h"
00012 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "FWCore/ParameterSet/interface/InputTag.h"
00017 #include "DataFormats/Candidate/interface/OverlapChecker.h"
00018 #include "TH1.h"
00019 #include <vector>
00020 
00021 class ZMuMuEfficiency : public edm::EDAnalyzer {
00022 public:
00023   ZMuMuEfficiency(const edm::ParameterSet& pset);
00024 private:
00025   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00026   virtual void endJob();
00027 
00028   edm::InputTag zMuTrack_, zMuTrackMatchMap_; 
00029   edm::InputTag zMuStandAlone_, zMuStandAloneMatchMap_;
00030   edm::InputTag muons_, muonMatchMap_, muonIso_;
00031   edm::InputTag tracks_, trackIso_;
00032   edm::InputTag standAlone_, standAloneIso_;
00033   double zMassMin_, zMassMax_, ptmin_, etamax_, isomax_;
00034   size_t nbinsPt_, nbinsEta_;
00035   reco::CandidateRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
00036   OverlapChecker overlap_;
00037 
00038   //histograms for measuring tracker efficiency
00039   TH1D *h_etaStandAlone_, *h_etaMuonOverlappedToStandAlone_; 
00040   TH1D *h_ptStandAlone_, *h_ptMuonOverlappedToStandAlone_; 
00041 
00042   //histograms for measuring standalone efficiency
00043   TH1D *h_etaTrack_, *h_etaMuonOverlappedToTrack_;
00044   TH1D *h_ptTrack_, *h_ptMuonOverlappedToTrack_;
00045 
00046   int numberOfMatchedZMuSta_, numberOfMatchedZMuTrack_;
00047   int numberOfOverlappedStandAlone_, numberOfOverlappedTracks_;
00048 };
00049 
00050 #include "FWCore/ServiceRegistry/interface/Service.h"
00051 #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
00052 #include "DataFormats/Common/interface/Handle.h"
00053 #include "DataFormats/Candidate/interface/Particle.h"
00054 #include "DataFormats/Candidate/interface/Candidate.h"
00055 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00056 #include "DataFormats/Candidate/interface/CandAssociation.h"
00057 #include <iostream>
00058 #include <iterator>
00059 #include <cmath>
00060 using namespace std;
00061 using namespace reco;
00062 using namespace edm;
00063 
00064 typedef CandDoubleAssociations IsolationCollection;
00065 
00066 ZMuMuEfficiency::ZMuMuEfficiency(const ParameterSet& pset) : 
00067   zMuTrack_(pset.getParameter<InputTag>("zMuTrack")), 
00068   zMuTrackMatchMap_(pset.getParameter<InputTag>("zMuTrackMatchMap")), 
00069   zMuStandAlone_(pset.getParameter<InputTag>("zMuStandAlone")), 
00070   zMuStandAloneMatchMap_(pset.getParameter<InputTag>("zMuStandAloneMatchMap")), 
00071   muons_(pset.getParameter<InputTag>("muons")), 
00072   muonMatchMap_(pset.getParameter<InputTag>("muonMatchMap")), 
00073   muonIso_(pset.getParameter<InputTag>("muonIso")), 
00074   tracks_(pset.getParameter<InputTag>("tracks")), 
00075   trackIso_(pset.getParameter<InputTag>("trackIso")), 
00076   standAlone_(pset.getParameter<InputTag>("standAlone")), 
00077   standAloneIso_(pset.getParameter<InputTag>("standAloneIso")), 
00078   zMassMin_(pset.getUntrackedParameter<double>("zMassMin")), 
00079   zMassMax_(pset.getUntrackedParameter<double>("zMassMax")), 
00080   ptmin_(pset.getUntrackedParameter<double>("ptmin")), 
00081   etamax_(pset.getUntrackedParameter<double>("etamax")),  
00082   isomax_(pset.getUntrackedParameter<double>("isomax")), 
00083   nbinsPt_(pset.getUntrackedParameter<size_t>("nbinsPt")), 
00084   nbinsEta_(pset.getUntrackedParameter<size_t>("nbinsEta")) {
00085   Service<TFileService> fs;
00086   TFileDirectory trackEffDir = fs->mkdir("TrackEfficiency");
00087   h_etaStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonEta", 
00088                                             "StandAlone #eta for Z -> #mu + standalone", 
00089                                             nbinsEta_, -etamax_, etamax_);
00090   h_etaMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAloneEta", 
00091                                                             "Global muon overlapped to standAlone #eta for Z -> #mu + sa", 
00092                                                             nbinsEta_, -etamax_, etamax_);
00093   h_ptStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonPt", 
00094                                            "StandAlone p_{t} for Z -> #mu + standalone", 
00095                                            nbinsPt_, ptmin_, 200);
00096   h_ptMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAlonePt", 
00097                                                            "Global muon overlapped to standAlone p_{t} for Z -> #mu + sa", 
00098                                                            nbinsPt_, ptmin_, 200);
00099   
00100   
00101   TFileDirectory standaloneEffDir = fs->mkdir("StandaloneEfficiency");
00102   h_etaTrack_ = standaloneEffDir.make<TH1D>("TrackMuonEta", 
00103                                             "Track #eta for Z -> #mu + track", 
00104                                             nbinsEta_, -etamax_, etamax_);
00105   h_etaMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackEta", 
00106                                                             "Global muon overlapped to track #eta for Z -> #mu + tk", 
00107                                                             nbinsEta_, -etamax_, etamax_);
00108   h_ptTrack_ = standaloneEffDir.make<TH1D>("TrackMuonPt", 
00109                                            "Track p_{t} for Z -> #mu + track", 
00110                                            nbinsPt_, ptmin_, 200);
00111   h_ptMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackPt", 
00112                                                            "Global muon overlapped to track p_{t} for Z -> #mu + tk", 
00113                                                            nbinsPt_, ptmin_, 200);
00114   
00115   numberOfMatchedZMuSta_ = 0;
00116   numberOfMatchedZMuTrack_ = 0;
00117   numberOfOverlappedStandAlone_ = 0;
00118   numberOfOverlappedTracks_ = 0;
00119 }
00120 
00121 void ZMuMuEfficiency::analyze(const Event& event, const EventSetup& setup) {
00122   Handle<CandidateCollection> zMuTrack;  
00123   Handle<CandMatchMap> zMuTrackMatchMap; //Map of Z made by Mu + Track
00124   Handle<CandidateCollection> zMuStandAlone; 
00125   Handle<CandMatchMap> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone
00126   Handle<CandidateCollection> muons; //Collection of Muons
00127   Handle<CandMatchMap> muonMatchMap; 
00128   Handle<IsolationCollection> muonIso; 
00129   Handle<CandidateCollection> tracks; //Collection of Tracks
00130   Handle<IsolationCollection> trackIso; 
00131   Handle<CandidateCollection> standAlone; //Collection of StandAlone
00132   Handle<IsolationCollection> standAloneIso; 
00133   
00134   event.getByLabel(zMuTrack_, zMuTrack); 
00135   event.getByLabel(zMuStandAlone_, zMuStandAlone); 
00136   event.getByLabel(muons_, muons); 
00137   event.getByLabel(tracks_, tracks); 
00138   event.getByLabel(standAlone_, standAlone); 
00139   
00140   //TRACK
00141   if (zMuStandAlone->size() > 0) {
00142     event.getByLabel(zMuStandAloneMatchMap_, zMuStandAloneMatchMap); 
00143     event.getByLabel(muonIso_, muonIso); 
00144     event.getByLabel(standAloneIso_, standAloneIso); 
00145     event.getByLabel(muonMatchMap_, muonMatchMap); 
00146     for(size_t i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates
00147       const Candidate & zMuStaCand = (*zMuStandAlone)[i]; //the candidate
00148       CandidateRef zMuStaCandRef(zMuStandAlone,i);
00149       bool isMatched = false;
00150       CandMatchMap::const_iterator zMuStaMapIt = zMuStandAloneMatchMap->find(zMuStaCandRef);
00151       if(zMuStaMapIt != zMuStandAloneMatchMap->end()) isMatched = true;
00152       CandidateRef dau0 = zMuStaCand.daughter(0)->masterClone().castTo<CandidateRef>();
00153       CandidateRef dau1 = zMuStaCand.daughter(1)->masterClone().castTo<CandidateRef>();
00154       
00155       // Cuts
00156       if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) && 
00157          (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) && 
00158          (zMuStaCand.mass() > zMassMin_) && (zMuStaCand.mass() < zMassMax_) && 
00159          (isMatched)) {
00160         CandidateRef standAloneCandRef(standAlone,0);
00161         if(dau0.id() == standAloneCandRef.id()) {
00162           standAloneMuonCandRef_ = dau0;
00163           globalMuonCandRef_ = dau1;
00164         }
00165         if(dau1.id()== standAloneCandRef.id()) {
00166           standAloneMuonCandRef_ = dau1;
00167           globalMuonCandRef_ = dau0;
00168         }
00169         //The Z daughters are already matched!
00170         const double globalMuonIsolation = (*muonIso)[globalMuonCandRef_];
00171         const double standAloneMuonIsolation = (*standAloneIso)[standAloneMuonCandRef_];
00172         
00173         if((globalMuonIsolation < isomax_) && (standAloneMuonIsolation < isomax_)) {
00174           numberOfMatchedZMuSta_++;
00175           h_etaStandAlone_->Fill(standAloneMuonCandRef_->eta()); //Denominator eta for measuring track efficiency
00176           h_ptStandAlone_->Fill(standAloneMuonCandRef_->pt());   //Denominator pt for measuring track eff
00177           
00178           for(size_t j = 0; j < muons->size() ; ++j) {
00179             const Candidate & muCand = (*muons)[j]; 
00180             CandidateRef muCandRef(muons, j); 
00181             CandMatchMap::const_iterator muonMapIt = muonMatchMap->find(muCandRef); 
00182             if((muonMapIt != muonMatchMap->end()) && (overlap_(*standAloneMuonCandRef_, muCand))) { 
00183               h_etaMuonOverlappedToStandAlone_->Fill(standAloneMuonCandRef_->eta()); //Numerator eta
00184               h_ptMuonOverlappedToStandAlone_->Fill(standAloneMuonCandRef_->pt());   //Numerator pt
00185               numberOfOverlappedTracks_++;
00186             }
00187           }
00188         }
00189       }
00190     }
00191   } //end loop on Candidate
00192   
00193   //STANDALONE
00194   if (zMuTrack->size() > 0) {
00195     event.getByLabel(zMuTrackMatchMap_, zMuTrackMatchMap); 
00196     event.getByLabel(muonIso_, muonIso); 
00197     event.getByLabel(trackIso_, trackIso); 
00198     event.getByLabel(muonMatchMap_, muonMatchMap); 
00199     for(size_t i = 0; i < zMuTrack->size(); ++i) { //loop on candidates
00200       const Candidate & zMuTrkCand = (*zMuTrack)[i]; //the candidate
00201       CandidateRef zMuTrkCandRef(zMuTrack,i);
00202       bool isMatched = false;
00203       CandMatchMap::const_iterator zMuTrkMapIt = zMuTrackMatchMap->find(zMuTrkCandRef);
00204       if(zMuTrkMapIt != zMuTrackMatchMap->end()) isMatched = true;
00205       CandidateRef dau0 = zMuTrkCand.daughter(0)->masterClone().castTo<CandidateRef>();
00206       CandidateRef dau1 = zMuTrkCand.daughter(1)->masterClone().castTo<CandidateRef>();
00207 
00208       // Cuts
00209       if ((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) && 
00210           (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta())< etamax_) && 
00211           (zMuTrkCand.mass() > zMassMin_) && (zMuTrkCand.mass() < zMassMax_) && 
00212           (isMatched)) {
00213         CandidateRef trackCandRef(tracks,0);
00214         if(dau0.id() == trackCandRef.id()) {
00215           trackMuonCandRef_ = dau0;
00216           globalMuonCandRef_ = dau1;
00217         }
00218         if(dau1.id() == trackCandRef.id()) {
00219           trackMuonCandRef_ = dau1;
00220           globalMuonCandRef_ = dau0;
00221         }
00222         //The Z daughters are already matched!
00223         const double globalMuonIsolation = (*muonIso)[globalMuonCandRef_];
00224         const double trackMuonIsolation = (*trackIso)[trackMuonCandRef_];
00225         
00226         if((globalMuonIsolation < isomax_) && (trackMuonIsolation < isomax_)) {
00227           numberOfMatchedZMuTrack_++;
00228           h_etaTrack_->Fill(trackMuonCandRef_->eta()); //Denominator eta Sta
00229           h_ptTrack_->Fill(trackMuonCandRef_->pt());   //Denominator pt Sta
00230           
00231           for(size_t j = 0; j < muons->size() ; ++j) {
00232             const Candidate & muCand = (*muons)[j];
00233             CandidateRef muCandRef(muons, j); 
00234             CandMatchMap::const_iterator muonMapIt = muonMatchMap->find(muCandRef); 
00235             if((muonMapIt != muonMatchMap->end()) && (overlap_(*trackMuonCandRef_, muCand))) { 
00236               h_etaMuonOverlappedToTrack_->Fill(trackMuonCandRef_->eta()); //Numerator sta eta
00237               h_ptMuonOverlappedToTrack_->Fill(trackMuonCandRef_->pt());   //Numerator sta pt
00238               numberOfOverlappedStandAlone_++;
00239             }
00240           }
00241         }
00242       }
00243     }
00244   } //end loop on Candidate  
00245 } 
00246 
00247 
00248 void ZMuMuEfficiency::endJob() {
00249   double efficiencySTA =(double)numberOfOverlappedStandAlone_/(double)numberOfMatchedZMuTrack_;
00250   double errorEff_STA = sqrt( efficiencySTA*(1 - efficiencySTA)/(double)numberOfMatchedZMuTrack_);
00251 
00252   double efficiencyTRACK =(double)numberOfOverlappedTracks_/(double)numberOfMatchedZMuSta_;
00253   double errorEff_TRACK = sqrt( efficiencyTRACK*(1 - efficiencyTRACK)/(double)numberOfMatchedZMuSta_);
00254 
00255   cout << "------------------------------------   Efficiency   ----------------------------- " << endl;
00256   cout << "numberOfOverlappedStandAlone = " << numberOfOverlappedStandAlone_ << endl; 
00257   cout << "numberOfMatchedZMuTrack = " << numberOfMatchedZMuTrack_ << endl; 
00258   cout << "numberOfOverlappedTracks = " << numberOfOverlappedTracks_ << endl; 
00259   cout << "numberOfMatchedZMuSta = " << numberOfMatchedZMuSta_ << endl; 
00260   cout << "Efficiency StandAlone = " << efficiencySTA << " +/- " << errorEff_STA << endl;
00261   cout << "Efficiency Track      = " << efficiencyTRACK << " +/- " << errorEff_TRACK << endl;
00262 }
00263   
00264 #include "FWCore/Framework/interface/MakerMacros.h"
00265 
00266 DEFINE_FWK_MODULE(ZMuMuEfficiency);
00267
ZMuMuEfficiency.cc
