Reference Manual

Go to the documentation of this file.
00001 /* \class ZMuMu_MCanalyzer
00002  * 
00003  * author: Davide Piccolo
00004  *
00005  * ZMuMu MC analyzer:
00006  * check muon reco efficiencies from MC truth, 
00007  *
00008  */
00009 
00010 #include "DataFormats/Common/interface/AssociationVector.h"
00011 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00012 #include "DataFormats/Candidate/interface/CandMatchMap.h" 
00013 #include "FWCore/Framework/interface/EDAnalyzer.h"
00014 #include "DataFormats/Candidate/interface/Particle.h"
00015 #include "DataFormats/Candidate/interface/Candidate.h"
00016 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00017 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00018 #include "DataFormats/MuonReco/interface/Muon.h"
00019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/Framework/interface/EventSetup.h"
00022 #include "FWCore/Utilities/interface/InputTag.h"
00023 #include "DataFormats/Candidate/interface/OverlapChecker.h"
00024 #include "DataFormats/Math/interface/deltaR.h"
00025 #include "DataFormats/PatCandidates/interface/Muon.h" 
00026 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
00027 #include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
00028 #include "DataFormats/Common/interface/AssociationVector.h"
00029 #include "DataFormats/PatCandidates/interface/PATObject.h"
00030 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00031 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
00032 #include "DataFormats/PatCandidates/interface/Isolation.h"
00033 #include "DataFormats/Common/interface/ValueMap.h"
00034 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
00035 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
00036 
00037 
00038 
00039 
00040 #include "TH1.h"
00041 #include "TH2.h"
00042 #include "TH3.h"
00043 #include <vector>
00044 using namespace edm;
00045 using namespace std;
00046 using namespace reco;
00047 using namespace reco;
00048 using namespace isodeposit;
00049 
00050 
00051 class ZMuMu_MCanalyzer : public edm::EDAnalyzer {
00052 public:
00053   ZMuMu_MCanalyzer(const edm::ParameterSet& pset);
00054 private:
00055   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00056   bool check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
00057   float getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
00058   float getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
00059   float getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
00060   Particle::LorentzVector getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
00061   virtual void endJob();
00062 
00063   edm::InputTag zMuMu_, zMuMuMatchMap_; 
00064   edm::InputTag zMuStandAlone_, zMuStandAloneMatchMap_;
00065   edm::InputTag zMuTrack_, zMuTrackMatchMap_; 
00066   edm::InputTag muons_, muonMatchMap_, muonIso_;
00067   edm::InputTag tracks_, trackIso_;
00068   edm::InputTag genParticles_;
00069 
00070   bool bothMuons_;
00071 
00072   double etamin_, etamax_, ptmin_, massMin_, massMax_, isoMax_;
00073 
00074   double ptThreshold_, etEcalThreshold_, etHcalThreshold_ ,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_;
00075  
00076   bool relativeIsolation_;
00077   string hltPath_;
00078   reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
00079   OverlapChecker overlap_;
00080 
00081   // general histograms
00082   TH1D *h_trackProbe_eta, *h_trackProbe_pt, *h_staProbe_eta, *h_staProbe_pt, *h_ProbeOk_eta, *h_ProbeOk_pt;
00083 
00084   // global counters
00085   int nGlobalMuonsMatched_passed;    // total number of global muons MC matched and passing cuts (and triggered)
00086   int nGlobalMuonsMatched_passedIso;    // total number of global muons MC matched and passing cuts including Iso
00087   int n2GlobalMuonsMatched_passedIso;    // total number of Z->2 global muons MC matched and passing cuts including Iso
00088   int nStaMuonsMatched_passedIso;       // total number of sta only muons MC matched and passing cuts including Iso
00089   int nTracksMuonsMatched_passedIso;    // total number of tracks only muons MC matched and passing cuts including Iso
00090   int n2GlobalMuonsMatched_passedIso2Trg;    // total number of Z->2 global muons MC matched and passing cuts including Iso and both triggered
00091   int nMu0onlyTriggered;               // n. of events zMuMu with mu0 only triggered
00092   int nMu1onlyTriggered;               // n. of events zMuMu with mu1 only triggered
00093 
00094   int nZMuMu_matched;               // n. of events zMuMu macthed
00095   int nZMuSta_matched;              // n  of events zMuSta macthed 
00096   int nZMuTrk_matched;              // n. of events zMuTrk mathced
00097 };
00098 
00099 
00100 template<typename T>
00101 double isolation(const T * t, double ptThreshold, double etEcalThreshold, double etHcalThreshold , double dRVetoTrk, double dRTrk, double dREcal , double dRHcal,  double alpha, double beta, bool relativeIsolation) {
00102   // on 34X: 
00103   const pat::IsoDeposit * trkIso = t->isoDeposit(pat::TrackIso);
00104   //  const pat::IsoDeposit * trkIso = t->trackerIsoDeposit();
00105   // on 34X 
00106   const pat::IsoDeposit * ecalIso = t->isoDeposit(pat::EcalIso);
00107   //  const pat::IsoDeposit * ecalIso = t->ecalIsoDeposit();
00108   //    on 34X 
00109   const pat::IsoDeposit * hcalIso = t->isoDeposit(pat::HcalIso);   
00110   //    const pat::IsoDeposit * hcalIso = t->hcalIsoDeposit();
00111 
00112   Direction dir = Direction(t->eta(), t->phi());
00113    
00114   pat::IsoDeposit::AbsVetos vetosTrk;
00115   vetosTrk.push_back(new ConeVeto( dir, dRVetoTrk ));
00116   vetosTrk.push_back(new ThresholdVeto( ptThreshold ));
00117     
00118   pat::IsoDeposit::AbsVetos vetosEcal;
00119   vetosEcal.push_back(new ConeVeto( dir, 0.));
00120   vetosEcal.push_back(new ThresholdVeto( etEcalThreshold ));
00121     
00122   pat::IsoDeposit::AbsVetos vetosHcal;
00123   vetosHcal.push_back(new ConeVeto( dir, 0. ));
00124   vetosHcal.push_back(new ThresholdVeto( etHcalThreshold ));
00125 
00126   double isovalueTrk = (trkIso->sumWithin(dRTrk,vetosTrk));
00127   double isovalueEcal = (ecalIso->sumWithin(dREcal,vetosEcal));
00128   double isovalueHcal = (hcalIso->sumWithin(dRHcal,vetosHcal));
00129     
00130 
00131   double iso = alpha*( ((1+beta)/2*isovalueEcal) + ((1-beta)/2*isovalueHcal) ) + ((1-alpha)*isovalueTrk) ;
00132   if(relativeIsolation) iso /= t->pt();
00133   return iso;
00134 }
00135 
00136 
00137 double candidateIsolation( const reco::Candidate* c, double ptThreshold, double etEcalThreshold, double etHcalThreshold , double dRVetoTrk, double dRTrk, double dREcal , double dRHcal,  double alpha, double beta, bool relativeIsolation) {
00138   const pat::Muon * mu = dynamic_cast<const pat::Muon *>(&*c->masterClone());
00139   if(mu != 0) return isolation(mu, ptThreshold, etEcalThreshold, etHcalThreshold ,dRVetoTrk, dRTrk, dREcal , dRHcal,  alpha, beta, relativeIsolation);
00140   const pat::GenericParticle * trk = dynamic_cast<const pat::GenericParticle*>(&*c->masterClone());
00141   if(trk != 0) return isolation(trk,  ptThreshold, etEcalThreshold, etHcalThreshold ,dRVetoTrk, dRTrk, dREcal ,
00142                                 dRHcal,  alpha, beta, relativeIsolation);
00143   throw edm::Exception(edm::errors::InvalidReference) 
00144     << "Candidate daughter #0 is neither pat::Muons nor pat::GenericParticle\n";      
00145   return -1;
00146 }
00147 
00148 
00149 
00150 
00151 #include "FWCore/ServiceRegistry/interface/Service.h"
00152 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00153 #include "DataFormats/Common/interface/Handle.h"
00154 #include "DataFormats/Candidate/interface/Particle.h"
00155 #include "DataFormats/Candidate/interface/Candidate.h"
00156 #include "DataFormats/Common/interface/ValueMap.h"
00157 #include "DataFormats/Candidate/interface/CandAssociation.h"
00158 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00159 #include "DataFormats/Math/interface/LorentzVector.h"
00160 #include "DataFormats/TrackReco/interface/Track.h"
00161 #include <iostream>
00162 #include <iterator>
00163 #include <cmath>
00164 using namespace std;
00165 using namespace reco;
00166 using namespace edm;
00167 
00168 
00169 typedef edm::ValueMap<float> IsolationCollection;
00170 
00171 ZMuMu_MCanalyzer::ZMuMu_MCanalyzer(const ParameterSet& pset) : 
00172   zMuMu_(pset.getParameter<InputTag>("zMuMu")), 
00173   zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")), 
00174   zMuStandAlone_(pset.getParameter<InputTag>("zMuStandAlone")), 
00175   zMuStandAloneMatchMap_(pset.getParameter<InputTag>("zMuStandAloneMatchMap")), 
00176   zMuTrack_(pset.getParameter<InputTag>("zMuTrack")), 
00177   zMuTrackMatchMap_(pset.getParameter<InputTag>("zMuTrackMatchMap")), 
00178   muons_(pset.getParameter<InputTag>("muons")), 
00179   tracks_(pset.getParameter<InputTag>("tracks")), 
00180   genParticles_(pset.getParameter<InputTag>( "genParticles" ) ),
00181 
00182   bothMuons_(pset.getParameter<bool>("bothMuons")), 
00183   etamin_(pset.getUntrackedParameter<double>("etamin")),  
00184   etamax_(pset.getUntrackedParameter<double>("etamax")),  
00185   ptmin_(pset.getUntrackedParameter<double>("ptmin")), 
00186   massMin_(pset.getUntrackedParameter<double>("zMassMin")), 
00187   massMax_(pset.getUntrackedParameter<double>("zMassMax")), 
00188   isoMax_(pset.getUntrackedParameter<double>("isomax")),
00189   ptThreshold_(pset.getUntrackedParameter<double>("ptThreshold")),
00190   etEcalThreshold_(pset.getUntrackedParameter<double>("etEcalThreshold")),
00191   etHcalThreshold_(pset.getUntrackedParameter<double>("etHcalThreshold")),
00192   dRVetoTrk_(pset.getUntrackedParameter<double>("deltaRVetoTrk")),
00193   dRTrk_(pset.getUntrackedParameter<double>("deltaRTrk")),
00194   dREcal_(pset.getUntrackedParameter<double>("deltaREcal")),
00195   dRHcal_(pset.getUntrackedParameter<double>("deltaRHcal")),
00196   alpha_(pset.getUntrackedParameter<double>("alpha")),
00197   beta_(pset.getUntrackedParameter<double>("beta")),
00198   relativeIsolation_(pset.getUntrackedParameter<bool>("relativeIsolation")),
00199   hltPath_(pset.getUntrackedParameter<std::string >("hltPath")) { 
00200   Service<TFileService> fs;
00201 
00202   // binning of entries array (at moment defined by hand and not in cfg file)
00203   double etaRange[8] = {-2.5,-2.,-1.2,-0.8,0.8,1.2,2.,2.5};
00204   double ptRange[4] = {20.,40.,60.,100.};
00205 
00206   // general histograms
00207   h_trackProbe_eta = fs->make<TH1D>("trackProbeEta","Eta of tracks",7,etaRange);
00208   h_trackProbe_pt = fs->make<TH1D>("trackProbePt","Pt of tracks",3,ptRange);
00209   h_staProbe_eta = fs->make<TH1D>("standAloneProbeEta","Eta of standAlone",7,etaRange);
00210   h_staProbe_pt = fs->make<TH1D>("standAloneProbePt","Pt of standAlone",3,ptRange);
00211   h_ProbeOk_eta = fs->make<TH1D>("probeOkEta","Eta of probe Ok",7,etaRange);
00212   h_ProbeOk_pt = fs->make<TH1D>("probeOkPt","Pt of probe ok",3,ptRange);
00213 
00214   // clear global counters
00215   nGlobalMuonsMatched_passed = 0;
00216   nGlobalMuonsMatched_passedIso = 0;
00217   n2GlobalMuonsMatched_passedIso = 0;
00218   nStaMuonsMatched_passedIso = 0;
00219   nTracksMuonsMatched_passedIso = 0;
00220   n2GlobalMuonsMatched_passedIso2Trg = 0;
00221   nMu0onlyTriggered = 0;
00222   nMu1onlyTriggered = 0;
00223   nZMuMu_matched = 0;
00224   nZMuSta_matched = 0;
00225   nZMuTrk_matched = 0;
00226 }
00227 
00228 void ZMuMu_MCanalyzer::analyze(const Event& event, const EventSetup& setup) {
00229   Handle<CandidateView> zMuMu;  
00230   Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global 
00231   Handle<CandidateView> zMuStandAlone;  
00232   Handle<GenParticleMatch> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone
00233   Handle<CandidateView> zMuTrack;  
00234   Handle<GenParticleMatch> zMuTrackMatchMap; //Map of Z made by Mu + Track
00235   Handle<CandidateView> muons; //Collection of Muons
00236   Handle<CandidateView> tracks; //Collection of Tracks
00237 
00238   Handle<GenParticleCollection> genParticles;  // Collection of Generatd Particles
00239   
00240   event.getByLabel(zMuMu_, zMuMu); 
00241   event.getByLabel(zMuStandAlone_, zMuStandAlone); 
00242   event.getByLabel(zMuTrack_, zMuTrack); 
00243   event.getByLabel(genParticles_, genParticles);
00244   event.getByLabel(muons_, muons); 
00245   event.getByLabel(tracks_, tracks); 
00246 
00247   /*
00248   cout << "*********  zMuMu         size : " << zMuMu->size() << endl;
00249   cout << "*********  zMuStandAlone size : " << zMuStandAlone->size() << endl;
00250   cout << "*********  zMuTrack      size : " << zMuTrack->size() << endl;
00251   cout << "*********  muons         size : " << muons->size() << endl;      
00252   cout << "*********  tracks        size : " << tracks->size() << endl;
00253   */
00254   //      std::cout<<"Run-> "<<event.id().run()<<std::endl;
00255   //      std::cout<<"Event-> "<<event.id().event()<<std::endl; 
00256 
00257 
00258   bool zMuMu_found = false;
00259 
00260   // loop on ZMuMu
00261   if (zMuMu->size() > 0 ) {
00262     event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap); 
00263     for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
00264       const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
00265       CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
00266 
00267       const Candidate * lep0 = zMuMuCand.daughter( 0 );
00268       const Candidate * lep1 = zMuMuCand.daughter( 1 );
00269       const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
00270       //      double trkiso0 = muonDau0.trackIso();
00271       const pat::Muon & muonDau1 = dynamic_cast<const pat::Muon &>(*lep1->masterClone());
00272       //double trkiso1 = muonDau1.trackIso();
00273 
00274       double iso0 = candidateIsolation(lep0,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00275       double iso1 = candidateIsolation(lep1,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00276 
00277       double pt0 = zMuMuCand.daughter(0)->pt();
00278       double pt1 = zMuMuCand.daughter(1)->pt();
00279       double eta0 = zMuMuCand.daughter(0)->eta();
00280       double eta1 = zMuMuCand.daughter(1)->eta();
00281       double mass = zMuMuCand.mass();
00282 
00283       // HLT match
00284       const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
00285         muonDau0.triggerObjectMatchesByPath( hltPath_ );
00286       const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
00287         muonDau1.triggerObjectMatchesByPath( hltPath_ );
00288 
00289       bool trig0found = false;
00290       bool trig1found = false;
00291       if( mu0HLTMatches.size()>0 )
00292         trig0found = true;
00293       if( mu1HLTMatches.size()>0 )
00294         trig1found = true;
00295 
00296       GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef];
00297       if(zMuMuMatch.isNonnull()) {  // ZMuMu matched
00298         zMuMu_found = true;
00299         nZMuMu_matched++;       
00300         if (pt0>ptmin_ && pt1>ptmin_ && abs(eta0)>etamin_ && abs(eta1) >etamin_ && abs(eta0)<etamax_ && abs(eta1) <etamax_ && mass >massMin_ && mass < massMax_ && (trig0found || trig1found)) { // kinematic and trigger cuts passed
00301           nGlobalMuonsMatched_passed++; // first global Muon passed kine cuts 
00302           nGlobalMuonsMatched_passed++; // second global muon passsed kine cuts
00303           if (iso0<isoMax_) nGlobalMuonsMatched_passedIso++;       // first global muon passed the iso cut
00304           if (iso1<isoMax_) nGlobalMuonsMatched_passedIso++;       // second global muon passed the iso cut
00305           if (iso0<isoMax_ && iso1<isoMax_) {
00306             n2GlobalMuonsMatched_passedIso++;  // both muons passed iso cut
00307             if (trig0found && trig1found) n2GlobalMuonsMatched_passedIso2Trg++;  // both muons have HLT
00308             if (trig0found && !trig1found) nMu0onlyTriggered++;
00309             if (trig1found && !trig0found) nMu1onlyTriggered++;
00310             // histograms vs eta and pt
00311             if (trig1found) {         // check efficiency of muon0 not imposing the trigger on it 
00312               h_trackProbe_eta->Fill(eta0);
00313               h_trackProbe_pt->Fill(pt0);
00314               h_staProbe_eta->Fill(eta0);
00315               h_staProbe_pt->Fill(pt0);
00316               h_ProbeOk_eta->Fill(eta0);
00317               h_ProbeOk_pt->Fill(pt0);
00318             }
00319             if (trig0found) {         // check efficiency of muon1 not imposing the trigger on it 
00320               h_trackProbe_eta->Fill(eta1);
00321               h_staProbe_eta->Fill(eta1);
00322               h_trackProbe_pt->Fill(pt1);
00323               h_staProbe_pt->Fill(pt1);
00324               h_ProbeOk_eta->Fill(eta1);
00325               h_ProbeOk_pt->Fill(pt1);
00326             }
00327           }
00328         }
00329       } // end MC match
00330 
00331     }  // end loop on ZMuMu cand
00332   }    // end if ZMuMu size > 0
00333 
00334   // loop on ZMuSta
00335   bool zMuSta_found = false;
00336   if (!zMuMu_found && zMuStandAlone->size() > 0 ) {
00337     event.getByLabel(zMuStandAloneMatchMap_, zMuStandAloneMatchMap); 
00338     for(unsigned int i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates
00339       const Candidate & zMuStandAloneCand = (*zMuStandAlone)[i]; //the candidate
00340       CandidateBaseRef zMuStandAloneCandRef = zMuStandAlone->refAt(i);
00341       GenParticleRef zMuStandAloneMatch = (*zMuStandAloneMatchMap)[zMuStandAloneCandRef];
00342 
00343       const Candidate * lep0 = zMuStandAloneCand.daughter( 0 );
00344       const Candidate * lep1 = zMuStandAloneCand.daughter( 1 );
00345       const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
00346       //double trkiso0 = muonDau0.trackIso();
00347       //      const pat::Muon & muonDau1 = dynamic_cast<const pat::Muon &>(*lep1->masterClone());
00348       //double trkiso1 = muonDau1.trackIso();
00349      
00350       double iso0 = candidateIsolation(lep0,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00351       double iso1 = candidateIsolation(lep1,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00352 
00353       double pt0 = zMuStandAloneCand.daughter(0)->pt();
00354       double pt1 = zMuStandAloneCand.daughter(1)->pt();
00355       double eta0 = zMuStandAloneCand.daughter(0)->eta();
00356       double eta1 = zMuStandAloneCand.daughter(1)->eta();
00357       double mass = zMuStandAloneCand.mass();
00358 
00359       // HLT match (check just dau0 the global)
00360       const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
00361         muonDau0.triggerObjectMatchesByPath( hltPath_ );
00362 
00363       bool trig0found = false;
00364       if( mu0HLTMatches.size()>0 )
00365         trig0found = true;
00366       
00367       if(zMuStandAloneMatch.isNonnull()) {  // ZMuStandAlone matched
00368         zMuSta_found = true;
00369         nZMuSta_matched++;      
00370         if (pt0>ptmin_ && pt1>ptmin_ && abs(eta0)>etamin_ && abs(eta1)>etamin_ && abs(eta0)<etamax_ && abs(eta1) <etamax_ && mass >massMin_ && 
00371             mass < massMax_ && iso0<isoMax_ && iso1 < isoMax_ && trig0found) { // all cuts and trigger passed
00372           nStaMuonsMatched_passedIso++;
00373           // histograms vs eta and pt
00374           h_staProbe_eta->Fill(eta1);
00375           h_staProbe_pt->Fill(pt1);
00376         }
00377       } // end MC match
00378     }  // end loop on ZMuStandAlone cand
00379   }    // end if ZMuStandAlone size > 0
00380 
00381 
00382   // loop on ZMuTrack
00383   bool zMuTrack_found = false;
00384   if (!zMuMu_found && !zMuSta_found && zMuTrack->size() > 0 ) {
00385     event.getByLabel(zMuTrackMatchMap_, zMuTrackMatchMap); 
00386     for(unsigned int i = 0; i < zMuTrack->size(); ++i) { //loop on candidates
00387       const Candidate & zMuTrackCand = (*zMuTrack)[i]; //the candidate
00388       CandidateBaseRef zMuTrackCandRef = zMuTrack->refAt(i);
00389       const Candidate * lep0 = zMuTrackCand.daughter( 0 );
00390       const Candidate * lep1 = zMuTrackCand.daughter( 1 );
00391       const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
00392       //double trkiso0 = muonDau0.trackIso();
00393       //const pat::GenericParticle & trackDau1 = dynamic_cast<const pat::GenericParticle &>(*lep1->masterClone());
00394       //double trkiso1 = trackDau1.trackIso();
00395       double iso0 = candidateIsolation(lep0,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00396       double iso1 = candidateIsolation(lep1,ptThreshold_, etEcalThreshold_, etHcalThreshold_,dRVetoTrk_, dRTrk_, dREcal_ , dRHcal_,  alpha_, beta_, relativeIsolation_);
00397 
00398 
00399       double pt0 = zMuTrackCand.daughter(0)->pt();
00400       double pt1 = zMuTrackCand.daughter(1)->pt();
00401       double eta0 = zMuTrackCand.daughter(0)->eta();
00402       double eta1 = zMuTrackCand.daughter(1)->eta();
00403       double mass = zMuTrackCand.mass();
00404 
00405       // HLT match (check just dau0 the global)
00406       const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
00407         muonDau0.triggerObjectMatchesByPath( hltPath_ );
00408 
00409       bool trig0found = false;
00410       if( mu0HLTMatches.size()>0 )
00411         trig0found = true;
00412 
00413       GenParticleRef zMuTrackMatch = (*zMuTrackMatchMap)[zMuTrackCandRef];
00414       if(zMuTrackMatch.isNonnull()) {  // ZMuTrack matched
00415         zMuTrack_found = true;
00416         nZMuTrk_matched++;
00417         if (pt0>ptmin_ && pt1>ptmin_ && abs(eta0)>etamin_ && abs(eta1)>etamin_ && abs(eta0)<etamax_ && abs(eta1) <etamax_ && mass >massMin_ && 
00418             mass < massMax_ && iso0<isoMax_ && iso1 < isoMax_ && trig0found) { // all cuts and trigger passed
00419           nTracksMuonsMatched_passedIso++;
00420           // histograms vs eta and pt
00421           h_trackProbe_eta->Fill(eta1);
00422           h_trackProbe_pt->Fill(pt1);
00423         }
00424       }  // end MC match
00425     }  // end loop on ZMuTrack cand
00426   }    // end if ZMuTrack size > 0
00427 
00428 }       // end analyze
00429 
00430 bool ZMuMu_MCanalyzer::check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
00431 {
00432   int partId0 = dauGen0->pdgId();
00433   int partId1 = dauGen1->pdgId();
00434   int partId2 = dauGen2->pdgId();
00435   bool muplusFound=false;
00436   bool muminusFound=false;
00437   bool ZFound=false;
00438   if (partId0==13 || partId1==13 || partId2==13) muminusFound=true;
00439   if (partId0==-13 || partId1==-13 || partId2==-13) muplusFound=true;
00440   if (partId0==23 || partId1==23 || partId2==23) ZFound=true;
00441   return muplusFound*muminusFound*ZFound;   
00442 }
00443  
00444 float ZMuMu_MCanalyzer::getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
00445 {
00446   int partId0 = dauGen0->pdgId();
00447   int partId1 = dauGen1->pdgId();
00448   int partId2 = dauGen2->pdgId();
00449   float ptpart=0.;
00450   if (partId0 == ipart) {
00451     for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
00452       const Candidate * dauMuGen = dauGen0->daughter(k);
00453       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00454         ptpart = dauMuGen->pt();
00455       }
00456     }
00457   }
00458   if (partId1 == ipart) {
00459     for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
00460       const Candidate * dauMuGen = dauGen1->daughter(k);
00461       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00462         ptpart = dauMuGen->pt();
00463       }
00464     }
00465   }
00466   if (partId2 == ipart) {
00467     for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
00468       const Candidate * dauMuGen = dauGen2->daughter(k);
00469       if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
00470         ptpart = dauMuGen->pt();
00471       }
00472     }
00473   }
00474   return ptpart;
00475 }
00476  
00477 float ZMuMu_MCanalyzer::getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
00478 {
00479   int partId0 = dauGen0->pdgId();
00480   int partId1 = dauGen1->pdgId();
00481   int partId2 = dauGen2->pdgId();
00482   float etapart=0.;
00483   if (partId0 == ipart) {
00484     for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
00485       const Candidate * dauMuGen = dauGen0->daughter(k);
00486       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00487         etapart = dauMuGen->eta();
00488       }
00489     }
00490   }
00491   if (partId1 == ipart) {
00492     for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
00493       const Candidate * dauMuGen = dauGen1->daughter(k);
00494       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00495         etapart = dauMuGen->eta();
00496       }
00497     }
00498   }
00499   if (partId2 == ipart) {
00500     for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
00501       const Candidate * dauMuGen = dauGen2->daughter(k);
00502       if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
00503         etapart = dauMuGen->eta();
00504       }
00505     }
00506   }
00507   return etapart;
00508 }
00509 
00510 float ZMuMu_MCanalyzer::getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
00511 {
00512   int partId0 = dauGen0->pdgId();
00513   int partId1 = dauGen1->pdgId();
00514   int partId2 = dauGen2->pdgId();
00515   float phipart=0.;
00516   if (partId0 == ipart) {
00517     for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
00518       const Candidate * dauMuGen = dauGen0->daughter(k);
00519       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00520         phipart = dauMuGen->phi();
00521       }
00522     }
00523   }
00524   if (partId1 == ipart) {
00525     for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
00526       const Candidate * dauMuGen = dauGen1->daughter(k);
00527       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00528         phipart = dauMuGen->phi();
00529       }
00530     }
00531   }
00532   if (partId2 == ipart) {
00533     for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
00534       const Candidate * dauMuGen = dauGen2->daughter(k);
00535       if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
00536         phipart = dauMuGen->phi();
00537       }
00538     }
00539   }
00540   return phipart;
00541 }
00542 
00543 Particle::LorentzVector ZMuMu_MCanalyzer::getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
00544 {
00545   int partId0 = dauGen0->pdgId();
00546   int partId1 = dauGen1->pdgId();
00547   int partId2 = dauGen2->pdgId();
00548   Particle::LorentzVector p4part(0.,0.,0.,0.);
00549   if (partId0 == ipart) {
00550     for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
00551       const Candidate * dauMuGen = dauGen0->daughter(k);
00552       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00553         p4part = dauMuGen->p4();
00554       }
00555     }
00556   }
00557   if (partId1 == ipart) {
00558     for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
00559       const Candidate * dauMuGen = dauGen1->daughter(k);
00560       if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
00561         p4part = dauMuGen->p4();
00562       }
00563     }
00564   }
00565   if (partId2 == ipart) {
00566     for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
00567       const Candidate * dauMuGen = dauGen2->daughter(k);
00568       if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
00569         p4part = dauMuGen->p4();
00570       }
00571     }
00572   }
00573   return p4part;
00574 }
00575  
00576 
00577 
00578 void ZMuMu_MCanalyzer::endJob() {
00579   
00580   
00581   double eff_Iso = double(nGlobalMuonsMatched_passedIso)/nGlobalMuonsMatched_passed;
00582   double err_effIso = sqrt(eff_Iso*(1-eff_Iso)/nGlobalMuonsMatched_passed);
00583  
00584   double n1_afterIso = 2*n2GlobalMuonsMatched_passedIso2Trg+nMu0onlyTriggered+nMu1onlyTriggered+nTracksMuonsMatched_passedIso;
00585   double n2_afterIso = 2*n2GlobalMuonsMatched_passedIso2Trg+nMu0onlyTriggered+nMu1onlyTriggered+nStaMuonsMatched_passedIso;
00586   double nGLB_afterIso = 2*n2GlobalMuonsMatched_passedIso2Trg+nMu0onlyTriggered+nMu1onlyTriggered;
00587   double effSta_afterIso = (2*n2GlobalMuonsMatched_passedIso2Trg+nMu0onlyTriggered+nMu1onlyTriggered)/n1_afterIso;
00588   double effTrk_afterIso = (2*n2GlobalMuonsMatched_passedIso2Trg+nMu0onlyTriggered+nMu1onlyTriggered)/n2_afterIso;
00589   double effHLT_afterIso = (2.* n2GlobalMuonsMatched_passedIso2Trg)/(2.* n2GlobalMuonsMatched_passedIso2Trg + nMu0onlyTriggered + nMu1onlyTriggered);
00590   double err_effHLT_afterIso= sqrt( effHLT_afterIso * (1 - effHLT_afterIso)/nGLB_afterIso);
00591   double err_effsta_afterIso = sqrt(effSta_afterIso*(1-effSta_afterIso)/n1_afterIso);
00592   double err_efftrk_afterIso = sqrt(effTrk_afterIso*(1-effTrk_afterIso)/n2_afterIso);
00593  
00594 
00595   cout << "------------------------------------  Counters  --------------------------------" << endl;
00596 
00597   cout << "number of events zMuMu matched " << nZMuMu_matched << endl;
00598   cout << "number of events zMuSta matched " << nZMuSta_matched << endl;
00599   cout << "number of events zMuTk matched " << nZMuTrk_matched << endl;
00600   cout << "number of events zMuMu with mu0 only triggered " << nMu0onlyTriggered << endl;
00601   cout << "number of events zMuMu with mu1 only triggered " << nMu1onlyTriggered << endl;
00602   cout << "=========================================" << endl;
00603   cout << "n. of global muons MC matched and passing cuts:           " << nGlobalMuonsMatched_passed << endl;
00604   cout << "n. of global muons MC matched and passing also Iso cut:       " << nGlobalMuonsMatched_passedIso << endl;
00605   cout << "n. of Z -> 2 global muons MC matched and passing ALL cuts:    " << n2GlobalMuonsMatched_passedIso << endl;
00606   cout << "n. of ZMuSta MC matched and passing ALL cuts:    " << nStaMuonsMatched_passedIso << endl;
00607   cout << "n. of ZmuTrck MC matched and passing ALL cuts:   " << nTracksMuonsMatched_passedIso << endl;
00608   cout << "n. of Z -> 2 global muons MC matched and passing ALL cuts and both triggered: " << n2GlobalMuonsMatched_passedIso2Trg << endl;
00609   cout << "=================================================================================" << endl;
00610   cout << "Iso efficiency: " << eff_Iso << " +/- " << err_effIso << endl;
00611   cout << "HLT efficiency: " << effHLT_afterIso << " +/- " << err_effHLT_afterIso << endl;
00612   cout << "eff StandAlone (after Isocut) : " << effSta_afterIso << "+/-" << err_effsta_afterIso << endl;
00613   cout << "eff Tracker (after Isocut)    : " << effTrk_afterIso << "+/-" << err_efftrk_afterIso << endl;
00614 
00615 }
00616   
00617 #include "FWCore/Framework/interface/MakerMacros.h"
00618 
00619 DEFINE_FWK_MODULE(ZMuMu_MCanalyzer);
00620
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMu_MCanalyzer.cc
