Reference Manual

Go to the documentation of this file.
00001 #include "DQM/Physics/src/HiggsDQM.h"
00002 
00003 #include <memory>
00004 
00005 // Framework
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "FWCore/Framework/interface/LuminosityBlock.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "FWCore/ParameterSet/interface/FileInPath.h"
00011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00012 #include "FWCore/Framework/interface/Run.h"
00013 #include "DataFormats/Provenance/interface/EventID.h"
00014 
00015 // Candidate handling
00016 #include "DataFormats/Candidate/interface/Candidate.h"
00017 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00018 #include "DataFormats/Candidate/interface/OverlapChecker.h"
00019 #include "DataFormats/Candidate/interface/CompositeCandidate.h"
00020 #include "DataFormats/Candidate/interface/CompositeCandidateFwd.h"
00021 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00022 #include <DataFormats/EgammaCandidates/interface/GsfElectron.h>
00023 #include <DataFormats/EgammaCandidates/interface/GsfElectronFwd.h>
00024 #include "DataFormats/MuonReco/interface/Muon.h"
00025 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00026 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00027 
00028 #include "DQMServices/Core/interface/DQMStore.h"
00029 #include "DQMServices/Core/interface/MonitorElement.h"
00030 
00031 //#include "HiggsAnalysis/HiggsToZZ4Leptons/plugins/HZZ4LeptonsElectronAssociationMap.h"
00032 //#include "HiggsAnalysis/HiggsToZZ4Leptons/plugins/HZZ4LeptonsMuonAssociationMap.h"
00033 
00034 // vertexing
00035 // Transient tracks
00036 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00037 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00038 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00039 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00041 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00042 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00043 
00044 // Vertex utilities
00045 #include "DataFormats/VertexReco/interface/Vertex.h"
00046 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00047 
00048 // Geometry
00049 #include "DataFormats/Math/interface/Vector3D.h"
00050 #include "DataFormats/Math/interface/LorentzVector.h"
00051 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00052 #include "DataFormats/Common/interface/AssociationVector.h"
00053 #include "DataFormats/Math/interface/Point3D.h"
00054 #include "DataFormats/Math/interface/deltaR.h"
00055 #include "DataFormats/Math/interface/deltaPhi.h"
00056 
00057 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00058 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00059 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
00060 #include "DataFormats/DetId/interface/DetId.h"
00061 #include "DataFormats/Common/interface/RefToBase.h"
00062 
00063 
00064 // RECO
00065 // Isolation
00066 //#include "HiggsAnalysis/HiggsToZZ4Leptons/plugins/CandidateHadIsolation.h"
00067 //#include "HiggsAnalysis/HiggsToZZ4Leptons/plugins/CandidateTkIsolation.h"
00068 #include "DataFormats/MuonReco/interface/Muon.h"
00069 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00070 #include "DataFormats/MuonReco/interface/MuonIsolation.h" 
00071 #include "DataFormats/TrackReco/interface/Track.h"
00072 #include <DataFormats/EgammaCandidates/interface/GsfElectron.h>
00073 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00074 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00075 #include "DataFormats/VertexReco/interface/Vertex.h"
00076 #include "DataFormats/JetReco/interface/CaloJet.h"
00077 
00078 //MET
00079 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00080 #include "DataFormats/METReco/interface/PFMETCollection.h"
00081 #include "DataFormats/METReco/interface/GenMETCollection.h"
00082 #include "DataFormats/METReco/interface/CaloMET.h"
00083 #include "DataFormats/METReco/interface/PFMET.h"
00084 #include "DataFormats/METReco/interface/MET.h"
00085 #include "DataFormats/METReco/interface/METCollection.h"
00086 
00087 #include "DataFormats/Math/interface/LorentzVector.h"
00088 #include "TLorentzVector.h"
00089 
00090 #include <iostream>
00091 #include <iomanip>
00092 #include <stdio.h>
00093 #include <string>
00094 #include <sstream>
00095 #include <math.h>
00096 
00097 using namespace edm;
00098 using namespace std;
00099 using namespace reco;
00100 
00101 struct SortCandByDecreasingPt {
00102   bool operator()( const Candidate &c1, const Candidate &c2) const {
00103     return c1.pt() > c2.pt();
00104   }
00105 };
00106 
00107 
00108 double HiggsDQM::Distance( const reco::Candidate & c1, const reco::Candidate & c2 ) {
00109         return  deltaR(c1,c2);
00110 }
00111 
00112 double HiggsDQM::DistancePhi( const reco::Candidate & c1, const reco::Candidate & c2 ) {
00113         return  deltaPhi(c1.p4().phi(),c2.p4().phi());
00114 }
00115 
00116 // This always returns only a positive deltaPhi
00117 double HiggsDQM::calcDeltaPhi(double phi1, double phi2) {
00118   double deltaPhi = phi1 - phi2;
00119   if (deltaPhi < 0) deltaPhi = -deltaPhi;
00120   if (deltaPhi > 3.1415926) {
00121     deltaPhi = 2 * 3.1415926 - deltaPhi;
00122   }
00123   return deltaPhi;
00124 }
00125 
00126 //
00127 // -- Constructor
00128 //
00129 HiggsDQM::HiggsDQM(const edm::ParameterSet& ps){
00130  //cout<<"Entering  HiggsDQM::HiggsDQM: "<<endl;
00131  
00132   edm::LogInfo("HZZ4LeptonsDQM") <<  " Creating HZZ4LeptonsDQM " << "\n" ;
00133   bei_ = Service<DQMStore>().operator->();
00134   typedef std::vector<edm::InputTag> vtag;
00135   // Get parameters from configuration file
00136   theElecTriggerPathToPass    = ps.getParameter<string>("elecTriggerPathToPass");
00137   theMuonTriggerPathToPass    = ps.getParameter<string>("muonTriggerPathToPass");
00138   theTriggerResultsCollection = ps.getParameter<InputTag>("triggerResultsCollection");
00139   theMuonCollectionLabel      = ps.getParameter<InputTag>("muonCollection");
00140   theElectronCollectionLabel  = ps.getParameter<InputTag>("electronCollection");
00141   theCaloJetCollectionLabel   = ps.getParameter<InputTag>("caloJetCollection");
00142   theCaloMETCollectionLabel   = ps.getParameter<InputTag>("caloMETCollection");
00143   thePfMETCollectionLabel     = ps.getParameter<InputTag>("pfMETCollection");
00144   // just to initialize
00145   isValidHltConfig_ = false;
00146   // cuts:
00147   ptThrMu1_ = ps.getUntrackedParameter<double>("PtThrMu1");
00148   ptThrMu2_ = ps.getUntrackedParameter<double>("PtThrMu2");
00149  
00150 
00151   
00152  //cout<<"...leaving  HiggsDQM::HiggsDQM. "<<endl;
00153 }
00154 //
00155 // -- Destructor
00156 //
00157 HiggsDQM::~HiggsDQM(){
00158   //cout<<"Entering HiggsDQM::~HiggsDQM: "<<endl;
00159   
00160   edm::LogInfo("HiggsDQM") <<  " Deleting HiggsDQM " << "\n" ;
00161 
00162   //cout<<"...leaving HiggsDQM::~HiggsDQM. "<<endl;
00163 }
00164 //
00165 // -- Begin Job
00166 //
00167 void HiggsDQM::beginJob(){
00168   //cout<<"Entering HiggsDQM::beginJob: "<<endl;
00169 
00170   nLumiSecs_ = 0;
00171   nEvents_   = 0;
00172   bei_->setCurrentFolder("Physics/Higgs");
00173   bookHistos(bei_);
00174   pi = 3.14159265;
00175   
00176 
00177   //cout<<"...leaving HiggsDQM::beginJob. "<<endl;
00178 }
00179 //
00180 // -- Begin Run
00181 //
00182 void HiggsDQM::beginRun(Run const& run, edm::EventSetup const& eSetup) {
00183   edm::LogInfo ("HiggsDQM") <<"[HiggsDQM]: Begining of Run";
00184   // passed as parameter to HLTConfigProvider::init(), not yet used
00185   bool isConfigChanged = false;
00186   
00187   // isValidHltConfig_ used to short-circuit analyze() in case of problems
00188   const std::string hltProcessName( "HLT" );
00189   isValidHltConfig_ = hltConfigProvider_.init( run, eSetup, hltProcessName, isConfigChanged );
00190 
00191 }
00192 //
00193 // -- Begin  Luminosity Block
00194 //
00195 void HiggsDQM::beginLuminosityBlock(edm::LuminosityBlock const& lumiSeg, 
00196                                             edm::EventSetup const& context) {
00197   //cout<<"Entering HiggsDQM::beginLuminosityBlock: "<<endl;
00198   
00199   edm::LogInfo ("HiggsDQM") <<"[HiggsDQM]: Begin of LS transition";
00200 
00201   //cout<<"...leaving HiggsDQM::beginLuminosityBlock. "<<endl;
00202 }
00203 //
00204 //  -- Book histograms
00205 //
00206 void HiggsDQM::bookHistos(DQMStore* bei){
00207   bei->cd();
00208   bei->setCurrentFolder("Physics/Higgs");
00209   h_vertex_number = bei->book1D("h_vertex_number", "Number of event vertices in collection", 10,-0.5,   9.5 );
00210   h_vertex_chi2  = bei->book1D("h_vertex_chi2" , "Event Vertex #chi^{2}/n.d.o.f."          , 100, 0.0,   2.0 );
00211   h_vertex_numTrks = bei->book1D("h_vertex_numTrks", "Event Vertex, number of tracks"     , 100, -0.5,  99.5 );
00212   h_vertex_sumTrks = bei->book1D("h_vertex_sumTrks", "Event Vertex, sum of track pt"      , 100,  0.0, 100.0 );
00213   h_vertex_d0    = bei->book1D("h_vertex_d0"   , "Event Vertex d0"                        , 100,  -10.0,  10.0);
00214   h_jet_et       = bei->book1D("h_jet_et",       "Jet with highest E_{T} (from "+theCaloJetCollectionLabel.label()+");E_{T}(1^{st} jet) (GeV)",    20, 0., 200.0);
00215   h_jet2_et      = bei->book1D("h_jet2_et",      "Jet with 2^{nd} highest E_{T} (from "+theCaloJetCollectionLabel.label()+");E_{T}(2^{nd} jet) (GeV)",    20, 0., 200.0);
00216   h_jet_count    = bei->book1D("h_jet_count",    "Number of "+theCaloJetCollectionLabel.label()+" (E_{T} > 15 GeV);Number of Jets", 8, -0.5, 7.5);
00217   h_caloMet      = bei->book1D("h_caloMet",        "Calo Missing E_{T}; GeV"          , 20,  0.0 , 100);
00218   h_caloMet_phi  = bei->book1D("h_caloMet_phi",    "Calo Missing E_{T} #phi;#phi(MET)", 35, -3.5, 3.5 );
00219   h_pfMet        = bei->book1D("h_pfMet",        "Pf Missing E_{T}; GeV"          , 20,  0.0 , 100);
00220   h_pfMet_phi    = bei->book1D("h_pfMet_phi",    "Pf Missing E_{T} #phi;#phi(MET)", 35, -3.5, 3.5 );
00221   h_eMultiplicity = bei_->book1D("NElectrons","# of electrons per event",10,0.,10.);
00222   h_mMultiplicity = bei_->book1D("NMuons","# of muons per event",10,0.,10.);
00223   h_ePt = bei_->book1D("ElePt","Pt of electrons",50,0.,100.);
00224   h_eEta = bei_->book1D("EleEta","Eta of electrons",100,-5.,5.);
00225   h_ePhi = bei_->book1D("ElePhi","Phi of electrons",100,-3.5,3.5);
00226   h_mPt_GMTM = bei_->book1D("MuonPt_GMTM","Pt of global+tracker muons",50,0.,100.);
00227   h_mEta_GMTM = bei_->book1D("MuonEta_GMTM","Eta of global+tracker muons",60,-3.,3.);
00228   h_mPhi_GMTM = bei_->book1D("MuonPhi_GMTM","Phi of global+tracker muons",70,-3.5,3.5);
00229   h_mPt_GMPT = bei_->book1D("MuonPt_GMPT","Pt of global prompt-tight muons",50,0.,100.);
00230   h_mEta_GMPT = bei_->book1D("MuonEta_GMPT","Eta of global prompt-tight muons",60,-3.,3.);
00231   h_mPhi_GMPT = bei_->book1D("MuonPhi_GMPT","Phi of global prompt-tight muons",70,-3.5,3.5);
00232   h_mPt_GM = bei_->book1D("MuonPt_GM","Pt of global muons",50,0.,100.);
00233   h_mEta_GM = bei_->book1D("MuonEta_GM","Eta of global muons",60,-3.,3.);
00234   h_mPhi_GM = bei_->book1D("MuonPhi_GM","Phi of global muons",70,-3.5,3.5);
00235   h_mPt_TM = bei_->book1D("MuonPt_TM","Pt of tracker muons",50,0.,100.);
00236   h_mEta_TM = bei_->book1D("MuonEta_TM","Eta of tracker muons",60,-3.,3.);
00237   h_mPhi_TM = bei_->book1D("MuonPhi_TM","Phi of tracker muons",70,-3.5,3.5);
00238   h_mPt_STAM = bei_->book1D("MuonPt_STAM","Pt of STA muons",50,0.,100.);
00239   h_mEta_STAM = bei_->book1D("MuonEta_STAM","Eta of STA muons",60,-3.,3.);
00240   h_mPhi_STAM = bei_->book1D("MuonPhi_STAM","Phi of STA muons",70,-3.5,3.5);
00241   h_eCombIso = bei_->book1D("EleCombIso","CombIso of electrons",100,0.,10.);
00242   h_mCombIso = bei_->book1D("MuonCombIso","CombIso of muons",100,0.,10.);
00243   h_dimumass_GMGM = bei->book1D("DimuMass_GMGM","Invariant mass of GMGM pairs",100,0.,200.);
00244   h_dimumass_GMTM = bei->book1D("DimuMass_GMTM","Invariant mass of GMTM pairs",100,0.,200.);
00245   h_dimumass_TMTM = bei->book1D("DimuMass_TMTM","Invariant mass of TMTM pairs",100,0.,200.);
00246   h_dielemass = bei->book1D("DieleMass","Invariant mass of EE pairs",100,0.,200.);
00247   h_lepcounts = bei->book1D("LeptonCounts","LeptonCounts for multi lepton events",49,0.,49.);
00248   
00249   bei->cd();  
00250 
00251 }
00252 //
00253 //  -- Analyze 
00254 //
00255 void HiggsDQM::analyze(const edm::Event& e, const edm::EventSetup& eSetup){
00256   //cout<<"[HiggsDQM::analyze()] "<<endl;
00257 
00258 //-------------------------------
00259 //--- Trigger Info
00260 //-------------------------------
00261   // short-circuit if hlt problems
00262   if( ! isValidHltConfig_ ) return;
00263   // Did it pass certain HLT path?
00264   Handle<TriggerResults> HLTresults;
00265   e.getByLabel(theTriggerResultsCollection, HLTresults); 
00266   if ( !HLTresults.isValid() ) return;
00267   //unsigned int triggerIndex_elec = hltConfig.triggerIndex(theElecTriggerPathToPass);
00268   //unsigned int triggerIndex_muon = hltConfig.triggerIndex(theMuonTriggerPathToPass);
00269   bool passed_electron_HLT = true;
00270   bool passed_muon_HLT     = true;
00271   //if (triggerIndex_elec < HLTresults->size()) passed_electron_HLT = HLTresults->accept(triggerIndex_elec);
00272   //if (triggerIndex_muon < HLTresults->size()) passed_muon_HLT     = HLTresults->accept(triggerIndex_muon);
00273   //if ( !(passed_electron_HLT || passed_muon_HLT) ) return;
00274 
00275 //-------------------------------
00276 //--- Vertex Info
00277 //-------------------------------
00278   Handle<VertexCollection> vertexHandle;
00279   e.getByLabel("offlinePrimaryVertices", vertexHandle);
00280   if ( vertexHandle.isValid() ){
00281     VertexCollection vertexCollection = *(vertexHandle.product());
00282     int vertex_number     = vertexCollection.size();
00283     VertexCollection::const_iterator v = vertexCollection.begin();
00284     double vertex_chi2    = v->normalizedChi2(); //v->chi2();
00285     double vertex_d0      = sqrt(v->x()*v->x()+v->y()*v->y());
00286     //double vertex_ndof    = v->ndof();cout << "ndof="<<vertex_ndof<<endl;
00287     double vertex_numTrks = v->tracksSize();
00288     double vertex_sumTrks = 0.0;
00289     for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack!=v->tracks_end(); vertex_curTrack++) {
00290       vertex_sumTrks += (*vertex_curTrack)->pt();
00291     }
00292     h_vertex_number->Fill(vertex_number);
00293     h_vertex_chi2->Fill(vertex_chi2);
00294     h_vertex_d0  ->Fill(vertex_d0);
00295     h_vertex_numTrks->Fill(vertex_numTrks);
00296     h_vertex_sumTrks->Fill(vertex_sumTrks);
00297   }
00298   
00299 //-------------------------------
00300 //--- Electrons
00301 //-------------------------------
00302   float nEle=0;
00303   Handle<GsfElectronCollection> electronCollection;
00304   e.getByLabel(theElectronCollectionLabel, electronCollection);
00305   if ( electronCollection.isValid() ){
00306     int posEle=0,negEle=0;
00307     // If it passed electron HLT and the collection was found, find electrons near Z mass
00308     if( passed_electron_HLT ) {
00309       for (reco::GsfElectronCollection::const_iterator recoElectron=electronCollection->begin(); recoElectron!=electronCollection->end(); recoElectron++){
00310 //      cout << "Electron with pt= " <<  recoElectron->pt() << " and eta" << recoElectron->eta() << " p=" <<  recoElectron->p() << endl;
00311         h_ePt->Fill(recoElectron->pt());
00312         h_eEta->Fill(recoElectron->eta());
00313         h_ePhi->Fill(recoElectron->phi());
00314         if(recoElectron->charge()==1){
00315           posEle++;
00316         }else if(recoElectron->charge()==-1){
00317           negEle++;
00318         }
00319       // Require electron to pass some basic cuts
00320       //if ( recoElectron->et() < 20 || fabs(recoElectron->eta())>2.5 ) continue;
00321       // Tighter electron cuts
00322       //if ( recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 || 
00323       //     recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 || 
00324       //     recoElectron->sigmaIetaIeta() > 0.027 ) continue;
00325       } // end of loop over electrons
00326     } // end if passed HLT
00327     nEle = posEle+negEle; if(nEle>9.) nEle=9.;
00328     h_eMultiplicity->Fill(nEle);  
00329 
00330     // Z->ee:
00331     unsigned int eleCollectionSize = electronCollection->size();
00332     for(unsigned int i=0; i<eleCollectionSize; i++) {
00333       const GsfElectron& ele = electronCollection->at(i);
00334       double pt = ele.pt();
00335       if(pt>ptThrMu1_){
00336         for(unsigned int j=i+1; j<eleCollectionSize; j++) {
00337           const GsfElectron& ele2 = electronCollection->at(j);
00338           double pt2 = ele2.pt();
00339           if(pt2>ptThrMu2_){
00340             const math::XYZTLorentzVector ZRecoEE (ele.px()+ele2.px(), ele.py()+ele2.py() , ele.pz()+ele2.pz(), ele.p()+ele2.p());
00341             h_dielemass->Fill(ZRecoEE.mass());
00342           }
00343         }
00344       }
00345     }
00346   }
00347   
00348   
00349   
00350 //-------------------------------
00351 //--- Muons
00352 //-------------------------------
00353   float nMu=0;
00354   Handle<MuonCollection> muonCollection;
00355   e.getByLabel(theMuonCollectionLabel,muonCollection);
00356   if ( muonCollection.isValid() ){
00357     // Find the highest pt muons
00358     int posMu=0,negMu=0;
00359     TLorentzVector m1, m2;
00360     if( passed_muon_HLT ) {
00361       for (reco::MuonCollection::const_iterator recoMuon=muonCollection->begin(); recoMuon!=muonCollection->end(); recoMuon++){
00362         //cout << "Muon with pt= " <<  muIter->pt() << " and eta" << muIter->eta() << " p=" <<  muIter->p() << endl;
00363         if(recoMuon->isGlobalMuon()&&recoMuon->isTrackerMuon()){
00364           h_mPt_GMTM->Fill(recoMuon->pt());
00365           h_mEta_GMTM->Fill(recoMuon->eta());
00366           h_mPhi_GMTM->Fill(recoMuon->phi());
00367         }else if(recoMuon->isGlobalMuon()&&(muon::isGoodMuon( (*recoMuon),muon::GlobalMuonPromptTight))){
00368           h_mPt_GMPT->Fill(recoMuon->pt());
00369           h_mEta_GMPT->Fill(recoMuon->eta());
00370           h_mPhi_GMPT->Fill(recoMuon->phi());
00371         }else if(recoMuon->isGlobalMuon()){
00372           h_mPt_GM->Fill(recoMuon->pt());
00373           h_mEta_GM->Fill(recoMuon->eta());
00374           h_mPhi_GM->Fill(recoMuon->phi());
00375         }else if(recoMuon->isTrackerMuon()&&(muon::segmentCompatibility((*recoMuon),reco::Muon::SegmentAndTrackArbitration))){
00376           h_mPt_TM->Fill(recoMuon->pt());
00377           h_mEta_TM->Fill(recoMuon->eta());
00378           h_mPhi_TM->Fill(recoMuon->phi());
00379         }else if(recoMuon->isStandAloneMuon()){
00380           h_mPt_STAM->Fill(recoMuon->pt());
00381           h_mEta_STAM->Fill(recoMuon->eta());
00382           h_mPhi_STAM->Fill(recoMuon->phi());
00383         }
00384         if ( recoMuon->charge()==1 ){
00385           posMu++;
00386         }else if ( recoMuon->charge()==-1 ){
00387           negMu++;
00388         }
00389       }
00390       nMu = posMu+negMu; if(nMu>9.) nMu=9.;
00391       h_mMultiplicity->Fill(nMu);
00392     }
00393 
00394     // Z->mumu:
00395     unsigned int muonCollectionSize = muonCollection->size();
00396     for(unsigned int i=0; i<muonCollectionSize; i++) {
00397       const Muon& mu = muonCollection->at(i);
00398       //if (!mu.isGlobalMuon()) continue;
00399       double pt = mu.pt();
00400       if(pt>ptThrMu1_){
00401         for(unsigned int j=i+1; j<muonCollectionSize; j++) {
00402           const Muon& mu2 = muonCollection->at(j);
00403           double pt2 = mu2.pt();
00404           if(pt2>ptThrMu2_){
00405             // Glb + Glb  
00406             if(mu.isGlobalMuon() && mu2.isGlobalMuon()){
00407               const math::XYZTLorentzVector ZRecoGMGM (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00408               h_dimumass_GMGM->Fill(ZRecoGMGM.mass());
00409             }
00410             // Glb + TM 
00411             else if(mu.isGlobalMuon() && mu2.isTrackerMuon()){
00412               const math::XYZTLorentzVector ZRecoGMTM (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00413               h_dimumass_GMTM->Fill(ZRecoGMTM.mass());
00414             }
00415             // TM + TM 
00416             else if(mu.isTrackerMuon() && mu2.isTrackerMuon()){
00417               const math::XYZTLorentzVector ZRecoTMTM (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00418               h_dimumass_TMTM->Fill(ZRecoTMTM.mass());
00419             }
00420           }
00421         }
00422       }
00423     }
00424   }
00425   
00426 //-------------------------------
00427 //--- Jets
00428 //-------------------------------
00429   Handle<CaloJetCollection> caloJetCollection;
00430   e.getByLabel (theCaloJetCollectionLabel,caloJetCollection);
00431   if ( caloJetCollection.isValid() ){
00432     float jet_et    = -8.0;
00433     float jet_eta   = -8.0;
00434     float jet_phi   = -8.0;
00435     int   jet_count = 0;
00436     float jet2_et   = -9.0;
00437     float jet2_eta  = -9.0;
00438     float jet2_phi  = -9.0;
00439     for (CaloJetCollection::const_iterator i_calojet = caloJetCollection->begin(); i_calojet != caloJetCollection->end(); i_calojet++) {
00440       float jet_current_et = i_calojet->et();
00441       // if it overlaps with electron, it is not a jet
00442       //if ( electron_et>0.0 && fabs(i_calojet->eta()-electron_eta ) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron_phi ) < 0.2) continue;
00443       //if ( electron2_et>0.0&& fabs(i_calojet->eta()-electron2_eta) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron2_phi) < 0.2) continue;
00444       // if it has too low Et, throw away
00445       if (jet_current_et < 15) continue;
00446       jet_count++;
00447       if (jet_current_et > jet_et) {
00448         jet2_et  = jet_et;  // 2nd highest jet get's et from current highest
00449         jet2_eta = jet_eta;
00450         jet2_phi = jet_phi;
00451         jet_et   = i_calojet->et(); // current highest jet gets et from the new highest
00452         jet_eta  = i_calojet->eta();
00453         jet_phi  = i_calojet->phi();
00454       } else if (jet_current_et > jet2_et) {
00455         jet2_et  = i_calojet->et();
00456         jet2_eta = i_calojet->eta();
00457         jet2_phi = i_calojet->phi();
00458       }
00459     }
00460     if (jet_et>0.0) {
00461       h_jet_et   ->Fill(jet_et);
00462       h_jet_count->Fill(jet_count);
00463     }
00464   }
00465   
00466 //-------------------------------
00467 //--- MET
00468 //-------------------------------
00469   Handle<CaloMETCollection> caloMETCollection;
00470   e.getByLabel(theCaloMETCollectionLabel, caloMETCollection);
00471   if ( caloMETCollection.isValid() ){
00472     float caloMet = caloMETCollection->begin()->et();
00473     float caloMet_phi = caloMETCollection->begin()->phi();
00474     h_caloMet        ->Fill(caloMet);
00475     h_caloMet_phi    ->Fill(caloMet_phi);
00476   }
00477   Handle<PFMETCollection> pfMETCollection;
00478   e.getByLabel(thePfMETCollectionLabel, pfMETCollection);
00479   if ( pfMETCollection.isValid() ){
00480     float pfMet = pfMETCollection->begin()->et();
00481     float pfMet_phi = pfMETCollection->begin()->phi();
00482     h_pfMet        ->Fill(pfMet);
00483     h_pfMet_phi    ->Fill(pfMet_phi);
00484   }
00485   
00486 //-------------------------------------
00487 //--- Events with more than 2 leptons:
00488 //-------------------------------------
00489   if(nMu+nEle > 2 && nMu+nEle < 10){
00490     if(nMu==0 && nEle==3) h_lepcounts->Fill(0);
00491     if(nMu==0 && nEle==4) h_lepcounts->Fill(1);
00492     if(nMu==0 && nEle==5) h_lepcounts->Fill(2);
00493     if(nMu==0 && nEle==6) h_lepcounts->Fill(3);
00494     if(nMu==0 && nEle==7) h_lepcounts->Fill(4);
00495     if(nMu==0 && nEle==8) h_lepcounts->Fill(5);
00496     if(nMu==0 && nEle==9) h_lepcounts->Fill(6);
00497     if(nMu==1 && nEle==2) h_lepcounts->Fill(7);
00498     if(nMu==1 && nEle==3) h_lepcounts->Fill(8);
00499     if(nMu==1 && nEle==4) h_lepcounts->Fill(9);
00500     if(nMu==1 && nEle==5) h_lepcounts->Fill(10);
00501     if(nMu==1 && nEle==6) h_lepcounts->Fill(11);
00502     if(nMu==1 && nEle==7) h_lepcounts->Fill(12);
00503     if(nMu==1 && nEle==8) h_lepcounts->Fill(13);
00504     if(nMu==2 && nEle==1) h_lepcounts->Fill(14);
00505     if(nMu==2 && nEle==2) h_lepcounts->Fill(15);
00506     if(nMu==2 && nEle==3) h_lepcounts->Fill(16);
00507     if(nMu==2 && nEle==4) h_lepcounts->Fill(17);
00508     if(nMu==2 && nEle==5) h_lepcounts->Fill(18);
00509     if(nMu==2 && nEle==6) h_lepcounts->Fill(19);
00510     if(nMu==2 && nEle==7) h_lepcounts->Fill(20);
00511     if(nMu==3 && nEle==0) h_lepcounts->Fill(21);
00512     if(nMu==3 && nEle==1) h_lepcounts->Fill(22);
00513     if(nMu==3 && nEle==2) h_lepcounts->Fill(23);
00514     if(nMu==3 && nEle==3) h_lepcounts->Fill(24);
00515     if(nMu==3 && nEle==4) h_lepcounts->Fill(25);
00516     if(nMu==3 && nEle==5) h_lepcounts->Fill(26);
00517     if(nMu==3 && nEle==6) h_lepcounts->Fill(27);
00518     if(nMu==4 && nEle==0) h_lepcounts->Fill(28);
00519     if(nMu==4 && nEle==1) h_lepcounts->Fill(29);
00520     if(nMu==4 && nEle==2) h_lepcounts->Fill(30);
00521     if(nMu==4 && nEle==3) h_lepcounts->Fill(31);
00522     if(nMu==4 && nEle==4) h_lepcounts->Fill(32);
00523     if(nMu==4 && nEle==5) h_lepcounts->Fill(33);
00524     if(nMu==5 && nEle==0) h_lepcounts->Fill(34);
00525     if(nMu==5 && nEle==1) h_lepcounts->Fill(35);
00526     if(nMu==5 && nEle==2) h_lepcounts->Fill(36);
00527     if(nMu==5 && nEle==3) h_lepcounts->Fill(37);
00528     if(nMu==5 && nEle==4) h_lepcounts->Fill(38);
00529     if(nMu==6 && nEle==0) h_lepcounts->Fill(39);
00530     if(nMu==6 && nEle==1) h_lepcounts->Fill(40);
00531     if(nMu==6 && nEle==2) h_lepcounts->Fill(41);
00532     if(nMu==6 && nEle==3) h_lepcounts->Fill(42);
00533     if(nMu==7 && nEle==0) h_lepcounts->Fill(43);
00534     if(nMu==7 && nEle==1) h_lepcounts->Fill(44);
00535     if(nMu==7 && nEle==2) h_lepcounts->Fill(45);
00536     if(nMu==8 && nEle==0) h_lepcounts->Fill(46);
00537     if(nMu==8 && nEle==1) h_lepcounts->Fill(47);
00538     if(nMu==9 && nEle==0) h_lepcounts->Fill(48);
00539     
00540   
00541   }
00542   if(nMu+nEle >= 10) cout<<"WARNING: "<<nMu+nEle<<" leptons in this event: run="<<e.id().run()<<", event="<<e.id().event()<<endl; 
00543   
00544 /*  /// channel conditions
00545   nElectron=0;
00546   nMuon=0;
00547   if (decaychannel=="2e2mu") {
00548     if ( posEle>=1 && negEle>=1 ) {
00549       nElectron=posEle+negEle;
00550     }
00551     if ( posMu>=1 && negMu>=1 ) {
00552       nMuon=posMu+negMu;
00553     }
00554   }
00555   else if (decaychannel=="4e") {
00556     if ( posEle>=2 && negEle>=2 ) {
00557       nElectron=posEle+negEle;
00558     }
00559   }
00560   else if (decaychannel=="4mu") {
00561     if ( posMu>=2 && negMu>=2 ) {
00562       nMuon=posMu+negMu;
00563     }  
00564   }
00565        
00566 
00567   // Pairs of EE MuMu 
00568   int nZEE=0,nZMuMu=0;
00569   if (decaychannel=="2e2mu"){
00570     Handle<CompositeCandidateCollection> zEECandidates;
00571     e.getByLabel(zToEETag_.label(), zEECandidates);    
00572     for ( CompositeCandidateCollection::const_iterator zIter=zEECandidates->begin(); zIter!= zEECandidates->end(); ++zIter ) {
00573 //      cout << "Zee mass= " << double(zIter->p4().mass()) << endl;
00574       if ( double(zIter->p4().mass())> 12. ){ 
00575         nZEE++;
00576       }
00577     }
00578     Handle<CompositeCandidateCollection> zMuMuCandidates;
00579     e.getByLabel(zToMuMuTag_.label(), zMuMuCandidates);
00580     for ( CompositeCandidateCollection::const_iterator zIter=zMuMuCandidates->begin(); zIter!= zMuMuCandidates->end(); ++zIter ) {
00581 //      cout << "Zmumu mass= " << double(zIter->p4().mass()) << endl;
00582       if ( zIter->p4().mass()> 12. ){
00583         nZMuMu++;
00584       }
00585     }    
00586   }
00587   
00588   // exclusive couples ZEE and ZMUMU
00589   if (decaychannel=="4e" ){
00590     Handle<CompositeCandidateCollection> higgsCandidates;
00591     e.getByLabel(hTozzTo4leptonsTag_.label(), higgsCandidates);
00592     for ( CompositeCandidateCollection::const_iterator hIter=higgsCandidates->begin(); hIter!= higgsCandidates->end(); ++hIter ) {
00593       if (nZEE<2) nZEE=0;
00594       for (size_t ndau=0; ndau<hIter->numberOfDaughters();ndau++){
00595         if ( hIter->daughter(ndau)->p4().mass()> 12.){  
00596           nZEE++;
00597         }
00598       }
00599     }
00600   }
00601   //
00602   if (decaychannel=="4mu" ){
00603     Handle<CompositeCandidateCollection> higgsCandidates;
00604     e.getByLabel(hTozzTo4leptonsTag_.label(), higgsCandidates);
00605     for ( CompositeCandidateCollection::const_iterator hIter=higgsCandidates->begin(); hIter!= higgsCandidates->end(); ++hIter ) {
00606       if (nZMuMu<2) nZMuMu=0;
00607       for (size_t ndau=0; ndau<hIter->numberOfDaughters();ndau++){
00608         if ( hIter->daughter(ndau)->p4().mass()> 12. ){  
00609           nZMuMu++;
00610         }
00611       }
00612     }
00613   }
00614 
00615   // 4 lepton combinations
00616   Handle<CompositeCandidateCollection> higgsCandidates;
00617   e.getByLabel(hTozzTo4leptonsTag_.label(), higgsCandidates);
00618    int nHiggs=0;
00619   for ( CompositeCandidateCollection::const_iterator hIter=higgsCandidates->begin(); hIter!= higgsCandidates->end(); ++hIter ) {
00620     if ( hIter->p4().mass()> 100. ){  
00621       nHiggs++;
00622     }
00623   }    
00624 */
00625 //  cout<<"Z and Higgs candidates: "<<nZEE<<" "<<nZMuMu<<" "<<nHiggs<<endl;
00626 
00627   //cout<<"[leaving HiggsDQM::analyze()] "<<endl;
00628 
00629 }
00630 //
00631 // -- End Luminosity Block
00632 //
00633 void HiggsDQM::endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, edm::EventSetup const& eSetup) {
00634 //  cout<<"Entering HiggsDQM::endLuminosityBlock: "<<endl;
00635 
00636   edm::LogInfo ("HiggsDQM") <<"[HiggsDQM]: End of LS transition, performing the DQM client operation";
00637 
00638   nLumiSecs_++;
00639   //cout << "nLumiSecs_: "<< nLumiSecs_ << endl;
00640   
00641   edm::LogInfo("HiggsDQM") << "====================================================== " << endl << " ===> Iteration # " << nLumiSecs_ << " " << lumiSeg.luminosityBlock() << endl  << "====================================================== " << endl;
00642 
00643 //  cout<<"...leaving HiggsDQM::endLuminosityBlock. "<<endl;
00644 }
00645 //
00646 // -- End Run
00647 //
00648 void HiggsDQM::endRun(edm::Run const& run, edm::EventSetup const& eSetup){
00649 //  cout<<"Entering HiggsDQM::endRun: "<<endl;
00650 
00651   //edm::LogVerbatim ("HiggsDQM") <<"[HiggsDQM]: End of Run, saving  DQM output ";
00652   //int iRun = run.run();
00653   
00654 //  cout<<"...leaving HiggsDQM::endRun. "<<endl;
00655 }
00656 
00657 //
00658 // -- End Job
00659 //
00660 void HiggsDQM::endJob(){
00661 //  cout<<"In HiggsDQM::endJob "<<endl;
00662   edm::LogInfo("HiggsDQM") <<"[HiggsDQM]: endjob called!";
00663 
00664 }
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/DQM/Physics/src/HiggsDQM.cc
