/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMuEfficiency.cc
C++ | 723 lines | 601 code | 76 blank | 46 comment | 171 complexity | e167282e0fc44917350a41f57757cd1f MD5 | raw file
1/* \class ZMuMuEfficiency 2 * 3 * author: Pasquale Noli 4 * revised by Salvatore di Guida 5 * revised for CSA08 by Davide Piccolo 6 * 7 * Efficiency of reconstruction tracker and muon Chamber 8 * 9 */ 10 11#include "DataFormats/Common/interface/AssociationVector.h" 12#include "DataFormats/HepMCCandidate/interface/GenParticle.h" 13#include "DataFormats/Candidate/interface/CandMatchMap.h" 14#include "FWCore/Framework/interface/EDAnalyzer.h" 15#include "DataFormats/Candidate/interface/Candidate.h" 16#include "DataFormats/Candidate/interface/CandidateFwd.h" 17#include "FWCore/ParameterSet/interface/ParameterSet.h" 18#include "FWCore/Framework/interface/Event.h" 19#include "FWCore/Framework/interface/EventSetup.h" 20#include "FWCore/Utilities/interface/InputTag.h" 21#include "DataFormats/Candidate/interface/OverlapChecker.h" 22#include "TH1.h" 23#include <vector> 24using namespace edm; 25using namespace std; 26using namespace reco; 27 28class ZMuMuEfficiency : public edm::EDAnalyzer { 29public: 30 ZMuMuEfficiency(const edm::ParameterSet& pset); 31private: 32 virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override; 33 bool check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 34 float getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 35 float getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 36 float getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 37 Particle::LorentzVector getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 38 virtual void endJob() override; 39 40 edm::InputTag zMuMu_, zMuMuMatchMap_; 41 edm::InputTag zMuTrack_, zMuTrackMatchMap_; 42 edm::InputTag zMuStandAlone_, zMuStandAloneMatchMap_; 43 edm::InputTag muons_, muonMatchMap_, muonIso_; 44 edm::InputTag tracks_, trackIso_; 45 edm::InputTag standAlone_, standAloneIso_; 46 edm::InputTag genParticles_; 47 48 double zMassMin_, zMassMax_, ptmin_, etamax_, isomax_; 49 unsigned int nbinsPt_, nbinsEta_; 50 reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_; 51 OverlapChecker overlap_; 52 53 //histograms for measuring tracker efficiency 54 TH1D *h_etaStandAlone_, *h_etaMuonOverlappedToStandAlone_; 55 TH1D *h_ptStandAlone_, *h_ptMuonOverlappedToStandAlone_; 56 57 //histograms for measuring standalone efficiency 58 TH1D *h_etaTrack_, *h_etaMuonOverlappedToTrack_; 59 TH1D *h_ptTrack_, *h_ptMuonOverlappedToTrack_; 60 61 //histograms for MC acceptance 62 TH1D *h_nZMCfound_; 63 TH1D *h_ZetaGen_, *h_ZptGen_, *h_ZmassGen_; 64 TH1D *h_muetaGen_, *h_muptGen_, *h_muIsoGen_; 65 TH1D *h_dimuonPtGen_, *h_dimuonMassGen_, *h_dimuonEtaGen_; 66 TH1D *h_ZetaGenPassed_, *h_ZptGenPassed_, *h_ZmassGenPassed_; 67 TH1D *h_muetaGenPassed_, *h_muptGenPassed_, *h_muIsoGenPassed_; 68 TH1D *h_dimuonPtGenPassed_, *h_dimuonMassGenPassed_, *h_dimuonEtaGenPassed_; 69 //histograms for invarian mass resolution 70 TH1D *h_DELTA_ZMuMuMassReco_dimuonMassGen_, *h_DELTA_ZMuStaMassReco_dimuonMassGen_, *h_DELTA_ZMuTrackMassReco_dimuonMassGen_; 71 72 int numberOfEventsWithZMuMufound, numberOfEventsWithZMuStafound; 73 int numberOfMatchedZMuSta_,numberOfMatchedSelectedZMuSta_; 74 int numberOfMatchedZMuMu_, numberOfMatchedSelectedZMuMu_; 75 int numberOfOverlappedStandAlone_, numberOfOverlappedTracks_, numberOfMatchedZMuTrack_notOverlapped; 76 int numberOfMatchedZMuTrack_exclusive ,numberOfMatchedSelectedZMuTrack_exclusive; 77 int numberOfMatchedZMuTrack_matchedZMuMu, numberOfMatchedZMuTrack_matchedSelectedZMuMu ; 78 int totalNumberOfevents, totalNumberOfZfound, totalNumberOfZPassed; 79 int noMCmatching, ZMuTrack_exclusive_1match, ZMuTrack_exclusive_morematch; 80 int ZMuTrackselected_exclusive_1match, ZMuTrackselected_exclusive_morematch; 81 int ZMuTrack_ZMuMu_1match, ZMuTrack_ZMuMu_2match, ZMuTrack_ZMuMu_morematch; 82 83 int n_zMuMufound_genZsele, n_zMuStafound_genZsele, n_zMuTrkfound_genZsele; 84}; 85 86#include "FWCore/ServiceRegistry/interface/Service.h" 87#include "CommonTools/UtilAlgos/interface/TFileService.h" 88#include "DataFormats/Common/interface/Handle.h" 89#include "DataFormats/Candidate/interface/Particle.h" 90#include "DataFormats/Candidate/interface/Candidate.h" 91#include "DataFormats/Common/interface/ValueMap.h" 92#include "DataFormats/Candidate/interface/CandAssociation.h" 93#include "DataFormats/HepMCCandidate/interface/GenParticle.h" 94#include <iostream> 95#include <iterator> 96#include <cmath> 97using namespace std; 98using namespace reco; 99using namespace edm; 100 101 102typedef edm::ValueMap<float> IsolationCollection; 103 104ZMuMuEfficiency::ZMuMuEfficiency(const ParameterSet& pset) : 105 zMuMu_(pset.getParameter<InputTag>("zMuMu")), 106 zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")), 107 zMuTrack_(pset.getParameter<InputTag>("zMuTrack")), 108 zMuTrackMatchMap_(pset.getParameter<InputTag>("zMuTrackMatchMap")), 109 zMuStandAlone_(pset.getParameter<InputTag>("zMuStandAlone")), 110 zMuStandAloneMatchMap_(pset.getParameter<InputTag>("zMuStandAloneMatchMap")), 111 muons_(pset.getParameter<InputTag>("muons")), 112 muonMatchMap_(pset.getParameter<InputTag>("muonMatchMap")), 113 muonIso_(pset.getParameter<InputTag>("muonIso")), 114 tracks_(pset.getParameter<InputTag>("tracks")), 115 trackIso_(pset.getParameter<InputTag>("trackIso")), 116 standAlone_(pset.getParameter<InputTag>("standAlone")), 117 standAloneIso_(pset.getParameter<InputTag>("standAloneIso")), 118 genParticles_(pset.getParameter<InputTag>( "genParticles" ) ), 119 120 zMassMin_(pset.getUntrackedParameter<double>("zMassMin")), 121 zMassMax_(pset.getUntrackedParameter<double>("zMassMax")), 122 ptmin_(pset.getUntrackedParameter<double>("ptmin")), 123 etamax_(pset.getUntrackedParameter<double>("etamax")), 124 isomax_(pset.getUntrackedParameter<double>("isomax")), 125 nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")), 126 nbinsEta_(pset.getUntrackedParameter<unsigned int>("nbinsEta")) { 127 Service<TFileService> fs; 128 TFileDirectory trackEffDir = fs->mkdir("TrackEfficiency"); 129 130 // tracker efficiency distributions 131 h_etaStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonEta", 132 "StandAlone #eta for Z -> #mu + standalone", 133 nbinsEta_, -etamax_, etamax_); 134 h_etaMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAloneEta", 135 "Global muon overlapped to standAlone #eta for Z -> #mu + sa", 136 nbinsEta_, -etamax_, etamax_); 137 h_ptStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonPt", 138 "StandAlone p_{t} for Z -> #mu + standalone", 139 nbinsPt_, ptmin_, 100); 140 h_ptMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAlonePt", 141 "Global muon overlapped to standAlone p_{t} for Z -> #mu + sa", 142 nbinsPt_, ptmin_, 100); 143 144 145 // StandAlone efficiency distributions 146 TFileDirectory standaloneEffDir = fs->mkdir("StandaloneEfficiency"); 147 h_etaTrack_ = standaloneEffDir.make<TH1D>("TrackMuonEta", 148 "Track #eta for Z -> #mu + track", 149 nbinsEta_, -etamax_, etamax_); 150 h_etaMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackEta", 151 "Global muon overlapped to track #eta for Z -> #mu + tk", 152 nbinsEta_, -etamax_, etamax_); 153 h_ptTrack_ = standaloneEffDir.make<TH1D>("TrackMuonPt", 154 "Track p_{t} for Z -> #mu + track", 155 nbinsPt_, ptmin_, 100); 156 h_ptMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackPt", 157 "Global muon overlapped to track p_{t} for Z -> #mu + tk", 158 nbinsPt_, ptmin_, 100); 159 160 161 // inv. mass resolution studies 162 TFileDirectory invMassResolutionDir = fs->mkdir("invriantMassResolution"); 163 h_DELTA_ZMuMuMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuMu_invMassResolution","zMuMu invariant Mass Resolution",50,-25,25); 164 h_DELTA_ZMuStaMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuSta_invMassResolution","zMuSta invariant Mass Resolution",50,-25,25); 165 h_DELTA_ZMuTrackMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuTrack_invMassResolution","zMuTrack invariant Mass Resolution",50,-25,25); 166 167 168 // generator level histograms 169 TFileDirectory genParticleDir = fs->mkdir("genParticle"); 170 h_nZMCfound_ = genParticleDir.make<TH1D>("NumberOfgeneratedZeta","n. of generated Z per event",4,-.5,3.5); 171 h_ZetaGen_ = genParticleDir.make<TH1D>("generatedZeta","#eta of generated Z",100,-5.,5.); 172 h_ZptGen_ = genParticleDir.make<TH1D>("generatedZpt","pt of generated Z",100,0.,200.); 173 h_ZmassGen_ = genParticleDir.make<TH1D>("generatedZmass","mass of generated Z",100,0.,200.); 174 h_muetaGen_ = genParticleDir.make<TH1D>("generatedMuonEta","#eta of generated muons from Z decay",100,-5.,5.); 175 h_muptGen_ = genParticleDir.make<TH1D>("generatedMuonpt","pt of generated muons from Z decay",100,0.,200.); 176 h_dimuonEtaGen_ = genParticleDir.make<TH1D>("generatedDimuonEta","#eta of generated dimuon",100,-5.,5.); 177 h_dimuonPtGen_ = genParticleDir.make<TH1D>("generatedDimuonPt","pt of generated dimuon",100,0.,200.); 178 h_dimuonMassGen_ = genParticleDir.make<TH1D>("generatedDimuonMass","mass of generated dimuon",100,0.,200.); 179 h_ZetaGenPassed_ = genParticleDir.make<TH1D>("generatedZeta_passed","#eta of generated Z after cuts",100,-5.,5.); 180 h_ZptGenPassed_ = genParticleDir.make<TH1D>("generatedZpt_passed","pt of generated Z after cuts",100,0.,200.); 181 h_ZmassGenPassed_ = genParticleDir.make<TH1D>("generatedZmass_passed","mass of generated Z after cuts",100,0.,200.); 182 h_muetaGenPassed_ = genParticleDir.make<TH1D>("generatedMuonEta_passed","#eta of generated muons from Z decay after cuts",100,-5.,5.); 183 h_muptGenPassed_ = genParticleDir.make<TH1D>("generatedMuonpt_passed","pt of generated muons from Z decay after cuts",100,0.,200.); 184 h_dimuonEtaGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonEta_passed","#eta of generated dimuon after cuts",100,-5.,5.); 185 h_dimuonPtGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonPt_passed","pt of generated dimuon after cuts",100,0.,200.); 186 h_dimuonMassGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonMass_passed","mass of generated dimuon after cuts",100,0.,200.); 187 // to insert isolation histograms .............. 188 189 numberOfEventsWithZMuMufound = 0; 190 numberOfEventsWithZMuStafound = 0; 191 numberOfMatchedZMuMu_ = 0; 192 numberOfMatchedSelectedZMuMu_ = 0; 193 numberOfMatchedZMuSta_ = 0; 194 numberOfMatchedSelectedZMuSta_ = 0; 195 numberOfMatchedZMuTrack_matchedZMuMu = 0; 196 numberOfMatchedZMuTrack_matchedSelectedZMuMu = 0; 197 numberOfMatchedZMuTrack_exclusive = 0; 198 numberOfMatchedSelectedZMuTrack_exclusive = 0; 199 numberOfOverlappedStandAlone_ = 0; 200 numberOfOverlappedTracks_ = 0; 201 numberOfMatchedZMuTrack_notOverlapped = 0; 202 noMCmatching = 0; 203 ZMuTrack_exclusive_1match = 0; 204 ZMuTrack_exclusive_morematch = 0; 205 ZMuTrackselected_exclusive_1match = 0; 206 ZMuTrackselected_exclusive_morematch = 0; 207 ZMuTrack_ZMuMu_1match = 0; 208 ZMuTrack_ZMuMu_2match = 0; 209 ZMuTrack_ZMuMu_morematch = 0; 210 211 n_zMuMufound_genZsele = 0; 212 n_zMuStafound_genZsele = 0; 213 n_zMuTrkfound_genZsele = 0; 214 215 // generator counters 216 totalNumberOfevents = 0; 217 totalNumberOfZfound = 0; 218 totalNumberOfZPassed = 0; 219 220} 221 222void ZMuMuEfficiency::analyze(const Event& event, const EventSetup& setup) { 223 Handle<CandidateView> zMuMu; 224 Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global 225 Handle<CandidateView> zMuTrack; 226 Handle<GenParticleMatch> zMuTrackMatchMap; //Map of Z made by Mu + Track 227 Handle<CandidateView> zMuStandAlone; 228 Handle<GenParticleMatch> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone 229 Handle<CandidateView> muons; //Collection of Muons 230 Handle<GenParticleMatch> muonMatchMap; 231 Handle<IsolationCollection> muonIso; 232 Handle<CandidateView> tracks; //Collection of Tracks 233 Handle<IsolationCollection> trackIso; 234 Handle<CandidateView> standAlone; //Collection of StandAlone 235 Handle<IsolationCollection> standAloneIso; 236 Handle<GenParticleCollection> genParticles; // Collection of Generatd Particles 237 238 event.getByLabel(zMuMu_, zMuMu); 239 event.getByLabel(zMuTrack_, zMuTrack); 240 event.getByLabel(zMuStandAlone_, zMuStandAlone); 241 event.getByLabel(muons_, muons); 242 event.getByLabel(tracks_, tracks); 243 event.getByLabel(standAlone_, standAlone); 244 event.getByLabel(genParticles_, genParticles); 245 246 cout << "********* zMuMu size : " << zMuMu->size() << endl; 247 cout << "********* zMuStandAlone size : " << zMuStandAlone->size() << endl; 248 cout << "********* zMuTrack size : " << zMuTrack->size() << endl; 249 cout << "********* muons size : " << muons->size()<< endl; 250 cout << "********* standAlone size : " << standAlone->size()<< endl; 251 cout << "********* tracks size : " << tracks->size()<< endl; 252 cout << "********* generated size : " << genParticles->size()<< endl; 253 254 255 // generator level distributions 256 257 int nZMCfound = 0; 258 totalNumberOfevents++; 259 int ngen = genParticles->size(); 260 bool ZMuMuMatchedfound = false; 261 bool ZMuMuMatchedSelectedfound = false; 262 bool ZMuStaMatchedfound = false; 263 //bool ZMuStaMatchedSelectedfound = false; 264 int ZMuTrackMatchedfound = 0; 265 int ZMuTrackMatchedSelected_exclusivefound = 0; 266 267 double dimuonMassGen = 0; 268 269 for (int i=0; i<ngen; i++) { 270 const Candidate &genCand = (*genParticles)[i]; 271 272 // if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0 273 // cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters() << " daughters" << endl; 274 if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0 275 if(genCand.numberOfDaughters() == 3) { // possible Z0 decays in mu+ mu-, the 3rd daughter is the same Z0 276 const Candidate * dauGen0 = genCand.daughter(0); 277 const Candidate * dauGen1 = genCand.daughter(1); 278 const Candidate * dauGen2 = genCand.daughter(2); 279 if (check_ifZmumu(dauGen0, dauGen1, dauGen2)) { // Z0 in mu+ mu- 280 totalNumberOfZfound++; 281 nZMCfound++; 282 bool checkpt = false; 283 bool checketa = false; 284 bool checkmass = false; 285 float mupluspt, muminuspt, mupluseta, muminuseta; 286 mupluspt = getParticlePt(-13,dauGen0,dauGen1,dauGen2); 287 muminuspt = getParticlePt(13,dauGen0,dauGen1,dauGen2); 288 mupluseta = getParticleEta(-13,dauGen0,dauGen1,dauGen2); 289 muminuseta = getParticleEta(13,dauGen0,dauGen1,dauGen2); 290 //float muplusphi = getParticlePhi(-13,dauGen0,dauGen1,dauGen2); 291 //float muminusphi = getParticlePhi(13,dauGen0,dauGen1,dauGen2); 292 Particle::LorentzVector pZ(0, 0, 0, 0); 293 Particle::LorentzVector muplusp4 = getParticleP4(-13,dauGen0,dauGen1,dauGen2); 294 Particle::LorentzVector muminusp4 = getParticleP4(13,dauGen0,dauGen1,dauGen2); 295 pZ = muplusp4 + muminusp4; 296 double dimuon_pt = sqrt(pZ.x()*pZ.x()+pZ.y()*pZ.y()); 297 double tan_theta_half = tan(atan(dimuon_pt/pZ.z())/2.); 298 double dimuon_eta = 0.; 299 if (tan_theta_half>0) dimuon_eta = -log(tan(tan_theta_half)); 300 if (tan_theta_half<=0) dimuon_eta = log(tan(-tan_theta_half)); 301 302 dimuonMassGen = pZ.mass(); // dimuon invariant Mass at Generator Level 303 304 h_ZmassGen_->Fill(genCand.mass()); 305 h_ZetaGen_->Fill(genCand.eta()); 306 h_ZptGen_->Fill(genCand.pt()); 307 h_dimuonMassGen_->Fill(pZ.mass()); 308 h_dimuonEtaGen_->Fill(dimuon_eta); 309 h_dimuonPtGen_->Fill(dimuon_pt); 310 h_muetaGen_->Fill(mupluseta); 311 h_muetaGen_->Fill(muminuseta); 312 h_muptGen_->Fill(mupluspt); 313 h_muptGen_->Fill(muminuspt); 314 // dimuon 4-momentum 315 // h_mDimuonMC->Fill(pZ.mass()); 316 // h_ZminusDimuonMassMC->Fill(genCand.mass()-pZ.mass()); 317 // h_DeltaPhiMC->Fill(deltaPhi(muplusphi,muminusphi)); 318 // if (dauGen2==23) float z_eta = dauGen2->eta(); 319 // if (dauGen2==23) float Zpt = dauGen2->pt(); 320 // h_DeltaPhivsZPtMC->Fill(DeltaPhi(muplusphi,muminusphi),ZPt); 321 322 if (mupluspt > ptmin_ && muminuspt > ptmin_) checkpt = true; 323 if (mupluseta < etamax_ && muminuseta < etamax_) checketa = true; 324 if (genCand.mass()>zMassMin_ && genCand.mass()<zMassMax_) checkmass = true; 325 if (checkpt && checketa && checkmass) { 326 totalNumberOfZPassed++; 327 h_ZmassGenPassed_->Fill(genCand.mass()); 328 h_ZetaGenPassed_->Fill(genCand.eta()); 329 h_ZptGenPassed_->Fill(genCand.pt()); 330 h_dimuonMassGenPassed_->Fill(pZ.mass()); 331 h_dimuonEtaGenPassed_->Fill(dimuon_eta); 332 h_dimuonPtGenPassed_->Fill(dimuon_pt); 333 h_muetaGenPassed_->Fill(mupluseta); 334 h_muetaGenPassed_->Fill(muminuseta); 335 h_muptGenPassed_->Fill(mupluspt); 336 h_muptGenPassed_->Fill(muminuspt); 337 338 if (zMuMu->size() > 0 ) { 339 n_zMuMufound_genZsele++; 340 } 341 else if (zMuStandAlone->size() > 0 ) { 342 n_zMuStafound_genZsele++; 343 } 344 else { 345 n_zMuTrkfound_genZsele++; 346 } 347 348 } 349 } 350 351 } 352 } 353 } 354 h_nZMCfound_->Fill(nZMCfound); // number of Z found in the event at generator level 355 356 //TRACK efficiency (conto numero di eventi Zmumu global e ZmuSta (ricorda che sono due campioni esclusivi) 357 358 if (zMuMu->size() > 0 ) { 359 numberOfEventsWithZMuMufound++; 360 event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap); 361 event.getByLabel(muonIso_, muonIso); 362 event.getByLabel(standAloneIso_, standAloneIso); 363 event.getByLabel(muonMatchMap_, muonMatchMap); 364 for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates 365 const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate 366 CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i); 367 bool isMatched = false; 368 GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef]; 369 370 if(zMuMuMatch.isNonnull()) { // ZMuMu matched 371 isMatched = true; 372 numberOfMatchedZMuMu_++; 373 } 374 CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone(); 375 CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone(); 376 if (isMatched) ZMuMuMatchedfound = true; 377 378 // Cuts 379 if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) && 380 (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) && 381 (zMuMuCand.mass() > zMassMin_) && (zMuMuCand.mass() < zMassMax_) && 382 (isMatched)) { 383 //The Z daughters are already matched! 384 const double globalMuonIsolation0 = (*muonIso)[dau0]; 385 const double globalMuonIsolation1 = (*muonIso)[dau1]; 386 if((globalMuonIsolation0 < isomax_) && (globalMuonIsolation1 < isomax_)) { // ZMuMu matched and selected by cuts 387 ZMuMuMatchedSelectedfound = true; 388 numberOfMatchedSelectedZMuMu_++; 389 h_etaStandAlone_->Fill(dau0->eta()); // StandAlone found dau0, eta 390 h_etaStandAlone_->Fill(dau1->eta()); // StandAlone found dau1, eta 391 h_etaMuonOverlappedToStandAlone_->Fill(dau0->eta()); // is global muon so dau0 is also found as a track, eta 392 h_etaMuonOverlappedToStandAlone_->Fill(dau1->eta()); // is global muon so dau1 is also found as a track, eta 393 h_ptStandAlone_->Fill(dau0->pt()); // StandAlone found dau0, pt 394 h_ptStandAlone_->Fill(dau1->pt()); // StandAlone found dau1, pt 395 h_ptMuonOverlappedToStandAlone_->Fill(dau0->pt()); // is global muon so dau0 is also found as a track, pt 396 h_ptMuonOverlappedToStandAlone_->Fill(dau1->pt()); // is global muon so dau1 is also found as a track, pt 397 398 h_etaTrack_->Fill(dau0->eta()); // Track found dau0, eta 399 h_etaTrack_->Fill(dau1->eta()); // Track found dau1, eta 400 h_etaMuonOverlappedToTrack_->Fill(dau0->eta()); // is global muon so dau0 is also found as a StandAlone, eta 401 h_etaMuonOverlappedToTrack_->Fill(dau1->eta()); // is global muon so dau1 is also found as a StandAlone, eta 402 h_ptTrack_->Fill(dau0->pt()); // Track found dau0, pt 403 h_ptTrack_->Fill(dau1->pt()); // Track found dau1, pt 404 h_ptMuonOverlappedToTrack_->Fill(dau0->pt()); // is global muon so dau0 is also found as a StandAlone, pt 405 h_ptMuonOverlappedToTrack_->Fill(dau1->pt()); // is global muon so dau1 is also found as a StandAlone, pt 406 407 h_DELTA_ZMuMuMassReco_dimuonMassGen_->Fill(zMuMuCand.mass()-dimuonMassGen); 408 // check that the two muons are matched . .per ora è solo un mio controllo 409 for(unsigned int j = 0; j < muons->size() ; ++j) { 410 CandidateBaseRef muCandRef = muons->refAt(j); 411 GenParticleRef muonMatch = (*muonMatchMap)[muCandRef]; 412 // if (muonMatch.isNonnull()) cout << "mu match n. " << j << endl; 413 } 414 } 415 } 416 } 417 } 418 419 if (zMuStandAlone->size() > 0) { 420 numberOfEventsWithZMuStafound++; 421 event.getByLabel(zMuStandAloneMatchMap_, zMuStandAloneMatchMap); 422 event.getByLabel(muonIso_, muonIso); 423 event.getByLabel(standAloneIso_, standAloneIso); 424 event.getByLabel(muonMatchMap_, muonMatchMap); 425 for(unsigned int i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates 426 const Candidate & zMuStaCand = (*zMuStandAlone)[i]; //the candidate 427 CandidateBaseRef zMuStaCandRef = zMuStandAlone->refAt(i); 428 bool isMatched = false; 429 GenParticleRef zMuStaMatch = (*zMuStandAloneMatchMap)[zMuStaCandRef]; 430 if(zMuStaMatch.isNonnull()) { // ZMuSta Macthed 431 isMatched = true; 432 ZMuStaMatchedfound = true; 433 numberOfMatchedZMuSta_++; 434 } 435 CandidateBaseRef dau0 = zMuStaCand.daughter(0)->masterClone(); 436 CandidateBaseRef dau1 = zMuStaCand.daughter(1)->masterClone(); 437 438 // Cuts 439 if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) && 440 (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) && 441 (zMuStaCand.mass() > zMassMin_) && (zMuStaCand.mass() < zMassMax_) && 442 (isMatched)) { 443 CandidateBaseRef standAloneMuonCandRef_, globalMuonCandRef_; 444 if(dau0->isGlobalMuon()) { 445 standAloneMuonCandRef_ = dau1; 446 globalMuonCandRef_ = dau0; 447 } 448 if(dau1->isGlobalMuon()) { 449 standAloneMuonCandRef_ = dau0; 450 globalMuonCandRef_ = dau1; 451 } 452 453 const double globalMuonIsolation = (*muonIso)[globalMuonCandRef_]; 454 const double standAloneMuonIsolation = (*standAloneIso)[standAloneMuonCandRef_]; 455 456 if((globalMuonIsolation < isomax_) && (standAloneMuonIsolation < isomax_)) { // ZMuSta matched and selected 457 //ZMuStaMatchedSelectedfound = true; 458 numberOfMatchedSelectedZMuSta_++; 459 h_etaStandAlone_->Fill(standAloneMuonCandRef_->eta()); //Denominator eta for measuring track efficiency 460 h_ptStandAlone_->Fill(standAloneMuonCandRef_->pt()); //Denominator pt for measuring track eff 461 h_DELTA_ZMuStaMassReco_dimuonMassGen_->Fill(zMuStaCand.mass()-dimuonMassGen); // differnce between ZMuSta reco and dimuon mass gen 462 463 } 464 } 465 } 466 } //end loop on Candidate 467 468 //STANDALONE efficiency 469 470 if (zMuTrack->size() > 0) { 471 event.getByLabel(zMuTrackMatchMap_, zMuTrackMatchMap); 472 event.getByLabel(muonIso_, muonIso); 473 event.getByLabel(trackIso_, trackIso); 474 event.getByLabel(muonMatchMap_, muonMatchMap); 475 for(unsigned int i = 0; i < zMuTrack->size(); ++i) { //loop on candidates 476 const Candidate & zMuTrkCand = (*zMuTrack)[i]; //the candidate 477 CandidateBaseRef zMuTrkCandRef = zMuTrack->refAt(i); 478 bool isMatched = false; 479 GenParticleRef zMuTrkMatch = (*zMuTrackMatchMap)[zMuTrkCandRef]; 480 if(zMuTrkMatch.isNonnull()) { 481 isMatched = true; 482 } 483 CandidateBaseRef dau0 = zMuTrkCand.daughter(0)->masterClone(); 484 CandidateBaseRef dau1 = zMuTrkCand.daughter(1)->masterClone(); 485 486 if (isMatched) { 487 ZMuTrackMatchedfound++; 488 if (ZMuMuMatchedfound) numberOfMatchedZMuTrack_matchedZMuMu++; 489 if (ZMuMuMatchedSelectedfound) numberOfMatchedZMuTrack_matchedSelectedZMuMu++; 490 if (!ZMuMuMatchedfound) numberOfMatchedZMuTrack_exclusive++; 491 } 492 // Cuts 493 if ((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) && 494 (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta())< etamax_) && 495 (zMuTrkCand.mass() > zMassMin_) && (zMuTrkCand.mass() < zMassMax_) && 496 (isMatched) && !ZMuMuMatchedfound && !ZMuStaMatchedfound ) { 497 498 // dau0 is always the global muon, dau1 is the track for ZMuTrack collection 499 const double globalMuonIsolation = (*muonIso)[dau0]; 500 const double trackMuonIsolation = (*trackIso)[dau1]; 501 if((globalMuonIsolation < isomax_) && (trackMuonIsolation < isomax_)) { // ZMuTRack matched - selected without ZMuMu found (exclusive) 502 numberOfMatchedSelectedZMuTrack_exclusive++; 503 ZMuTrackMatchedSelected_exclusivefound++; 504 h_etaTrack_->Fill(dau1->eta()); //Denominator eta Sta 505 h_ptTrack_->Fill(dau1->pt()); //Denominator pt Sta 506 h_DELTA_ZMuTrackMassReco_dimuonMassGen_->Fill(zMuTrkCand.mass()-dimuonMassGen); 507 } 508 509 } 510 } 511 } //end loop on Candidate 512 513 if (!ZMuMuMatchedfound && !ZMuStaMatchedfound && ZMuTrackMatchedfound == 0) noMCmatching++; 514 if (!ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_exclusive_1match++; 515 if (!ZMuMuMatchedfound && ZMuTrackMatchedfound > 1) ZMuTrack_exclusive_morematch++; 516 if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound == 1) ZMuTrackselected_exclusive_1match++; 517 if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound > 1) ZMuTrackselected_exclusive_morematch++; 518 if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_ZMuMu_1match++; 519 if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 2) ZMuTrack_ZMuMu_2match++; 520 if (ZMuMuMatchedfound && ZMuTrackMatchedfound > 2) ZMuTrack_ZMuMu_morematch++; 521 522} 523 524bool ZMuMuEfficiency::check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2) 525{ 526 int partId0 = dauGen0->pdgId(); 527 int partId1 = dauGen1->pdgId(); 528 int partId2 = dauGen2->pdgId(); 529 bool muplusFound=false; 530 bool muminusFound=false; 531 bool ZFound=false; 532 if (partId0==13 || partId1==13 || partId2==13) muminusFound=true; 533 if (partId0==-13 || partId1==-13 || partId2==-13) muplusFound=true; 534 if (partId0==23 || partId1==23 || partId2==23) ZFound=true; 535 return muplusFound*muminusFound*ZFound; 536} 537 538float ZMuMuEfficiency::getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2) 539{ 540 int partId0 = dauGen0->pdgId(); 541 int partId1 = dauGen1->pdgId(); 542 int partId2 = dauGen2->pdgId(); 543 float ptpart=0.; 544 if (partId0 == ipart) { 545 for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) { 546 const Candidate * dauMuGen = dauGen0->daughter(k); 547 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 548 ptpart = dauMuGen->pt(); 549 } 550 } 551 } 552 if (partId1 == ipart) { 553 for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) { 554 const Candidate * dauMuGen = dauGen1->daughter(k); 555 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 556 ptpart = dauMuGen->pt(); 557 } 558 } 559 } 560 if (partId2 == ipart) { 561 for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) { 562 const Candidate * dauMuGen = dauGen2->daughter(k); 563 if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) { 564 ptpart = dauMuGen->pt(); 565 } 566 } 567 } 568 return ptpart; 569} 570 571float ZMuMuEfficiency::getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2) 572{ 573 int partId0 = dauGen0->pdgId(); 574 int partId1 = dauGen1->pdgId(); 575 int partId2 = dauGen2->pdgId(); 576 float etapart=0.; 577 if (partId0 == ipart) { 578 for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) { 579 const Candidate * dauMuGen = dauGen0->daughter(k); 580 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 581 etapart = dauMuGen->eta(); 582 } 583 } 584 } 585 if (partId1 == ipart) { 586 for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) { 587 const Candidate * dauMuGen = dauGen1->daughter(k); 588 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 589 etapart = dauMuGen->eta(); 590 } 591 } 592 } 593 if (partId2 == ipart) { 594 for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) { 595 const Candidate * dauMuGen = dauGen2->daughter(k); 596 if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) { 597 etapart = dauMuGen->eta(); 598 } 599 } 600 } 601 return etapart; 602} 603 604float ZMuMuEfficiency::getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2) 605{ 606 int partId0 = dauGen0->pdgId(); 607 int partId1 = dauGen1->pdgId(); 608 int partId2 = dauGen2->pdgId(); 609 float phipart=0.; 610 if (partId0 == ipart) { 611 for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) { 612 const Candidate * dauMuGen = dauGen0->daughter(k); 613 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 614 phipart = dauMuGen->phi(); 615 } 616 } 617 } 618 if (partId1 == ipart) { 619 for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) { 620 const Candidate * dauMuGen = dauGen1->daughter(k); 621 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 622 phipart = dauMuGen->phi(); 623 } 624 } 625 } 626 if (partId2 == ipart) { 627 for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) { 628 const Candidate * dauMuGen = dauGen2->daughter(k); 629 if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) { 630 phipart = dauMuGen->phi(); 631 } 632 } 633 } 634 return phipart; 635} 636 637Particle::LorentzVector ZMuMuEfficiency::getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2) 638{ 639 int partId0 = dauGen0->pdgId(); 640 int partId1 = dauGen1->pdgId(); 641 int partId2 = dauGen2->pdgId(); 642 Particle::LorentzVector p4part(0.,0.,0.,0.); 643 if (partId0 == ipart) { 644 for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) { 645 const Candidate * dauMuGen = dauGen0->daughter(k); 646 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 647 p4part = dauMuGen->p4(); 648 } 649 } 650 } 651 if (partId1 == ipart) { 652 for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) { 653 const Candidate * dauMuGen = dauGen1->daughter(k); 654 if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) { 655 p4part = dauMuGen->p4(); 656 } 657 } 658 } 659 if (partId2 == ipart) { 660 for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) { 661 const Candidate * dauMuGen = dauGen2->daughter(k); 662 if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) { 663 p4part = dauMuGen->p4(); 664 } 665 } 666 } 667 return p4part; 668} 669 670 671 672void ZMuMuEfficiency::endJob() { 673 // double efficiencySTA =(double)numberOfOverlappedStandAlone_/(double)numberOfMatchedZMuTrack_; 674 // double errorEff_STA = sqrt( efficiencySTA*(1 - efficiencySTA)/(double)numberOfMatchedZMuTrack_); 675 676 double myTrackEff = 2.*numberOfMatchedSelectedZMuMu_/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_); 677 double myErrTrackEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_)); 678 679 double myStaEff = 2.*numberOfMatchedSelectedZMuMu_/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuTrack_exclusive); 680 double myErrStaEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuTrack_exclusive)); 681 682 // double efficiencyTRACK =(double)numberOfOverlappedTracks_/(double)numberOfMatchedZMuSta_; 683 // double errorEff_TRACK = sqrt( efficiencyTRACK*(1 - efficiencyTRACK)/(double)numberOfMatchedZMuSta_); 684 685 cout << "------------------------------------ Counters for MC acceptance --------------------------------" << endl; 686 cout << "totalNumberOfevents = " << totalNumberOfevents << endl; 687 cout << "totalNumberOfZfound = " << totalNumberOfZfound << endl; 688 cout << "totalNumberOfZpassed = " << totalNumberOfZPassed << endl; 689 cout << "n. of events zMuMu found (gen level selected)" << n_zMuMufound_genZsele << endl; 690 cout << "n. of events zMuSta found (gen level selected)" << n_zMuStafound_genZsele << endl; 691 cout << "n. of events zMuTrk found (gen level selected)" << n_zMuTrkfound_genZsele << endl; 692 693 cout << "---------------------------- Counter for MC truth efficiency calculation--------------------- " << endl; 694 695 cout << "number of events with ZMuMu found = " << numberOfEventsWithZMuMufound << endl; 696 cout << "number of events with ZMuSta found = " << numberOfEventsWithZMuStafound << endl; 697 cout << "-------------------------------------------------------------------------------------- " << endl; 698 699 cout << "number of events without MC maching = " << noMCmatching << endl; 700 cout << "number of ZMuTrack exclsive 1 match = " << ZMuTrack_exclusive_1match << endl; 701 cout << "number of ZMuTrack exclsive more match = " << ZMuTrack_exclusive_morematch << endl; 702 cout << "number of ZMuTrack selected exclusive 1 match = " << ZMuTrackselected_exclusive_1match << endl; 703 cout << "number of ZMuTrack selected exclusive more match = " << ZMuTrackselected_exclusive_morematch << endl; 704 cout << "number of ZMuTrack ZMuMu 1 match = " << ZMuTrack_ZMuMu_1match << endl; 705 cout << "number of ZMuTrack ZMuMu 2 match = " << ZMuTrack_ZMuMu_2match << endl; 706 cout << "number of ZMuTrack ZMuMu more match = " << ZMuTrack_ZMuMu_morematch << endl; 707 cout << "numberOfMatchedZMuMu = " << numberOfMatchedZMuMu_ << endl; 708 cout << "numberOfMatchedSelectdZMuMu = " << numberOfMatchedSelectedZMuMu_ << endl; 709 cout << "numberOfMatchedZMuSta = " << numberOfMatchedZMuSta_ << endl; 710 cout << "numberOfMatchedSelectedZMuSta = " << numberOfMatchedSelectedZMuSta_ << endl; 711 cout << "numberOfMatchedZMuTrack_matchedZMuMu = " << numberOfMatchedZMuTrack_matchedZMuMu << endl; 712 cout << "numberOfMatchedZMuTrack_matchedSelectedZMuMu = " << numberOfMatchedZMuTrack_matchedSelectedZMuMu << endl; 713 cout << "numberOfMatchedZMuTrack exclusive = " << numberOfMatchedZMuTrack_exclusive << endl; 714 cout << "numberOfMatchedSelectedZMuTrack exclusive = " << numberOfMatchedSelectedZMuTrack_exclusive << endl; 715 cout << " ----------------------------- Efficiency --------------------------------- " << endl; 716 cout << "Efficiency StandAlone = " << myStaEff << " +/- " << myErrStaEff << endl; 717 cout << "Efficiency Track = " << myTrackEff << " +/- " << myErrTrackEff << endl; 718} 719 720#include "FWCore/Framework/interface/MakerMacros.h" 721 722DEFINE_FWK_MODULE(ZMuMuEfficiency); 723