Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2011/02/18 13:57:46 $
00005  *  $Revision: 1.16 $
00006  *  \author Michael B. Anderson, University of Wisconsin-Madison
00007  *  \author Will Parker, University of Wisconsin-Madison
00008  */
00009 
00010 #include "DQM/Physics/src/EwkDQM.h"
00011 #include "DataFormats/Candidate/interface/Candidate.h"
00012 #include "FWCore/Framework/interface/Event.h"
00013 #include "FWCore/Framework/interface/EventSetup.h"
00014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00015 
00016 #include "DQMServices/Core/interface/DQMStore.h"
00017 #include "DQMServices/Core/interface/MonitorElement.h"
00018 #include "FWCore/ServiceRegistry/interface/Service.h"
00019 
00020 #include "DataFormats/Common/interface/Handle.h"
00021 
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 
00024 // Physics Objects
00025 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00026 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00027 #include "DataFormats/MuonReco/interface/Muon.h"
00028 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00029 #include "DataFormats/JetReco/interface/CaloJet.h"
00030 #include "DataFormats/METReco/interface/CaloMET.h"
00031 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00032 #include "DataFormats/METReco/interface/CaloMETFwd.h"
00033 #include "DataFormats/TrackReco/interface/Track.h"
00034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00035 #include "DataFormats/VertexReco/interface/Vertex.h"
00036 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00037 
00038 
00039 #include "DataFormats/Math/interface/LorentzVector.h"
00040 #include "TLorentzVector.h"
00041 
00042 #include <vector>
00043 
00044 #include <string>
00045 #include <cmath>
00046 using namespace std;
00047 using namespace edm;
00048 using namespace reco;
00049 
00050 
00051 
00052 EwkDQM::EwkDQM(const ParameterSet& parameters) {
00053   // Get parameters from configuration file
00054   theElecTriggerPathToPass    = parameters.getParameter<string>("elecTriggerPathToPass");
00055   theMuonTriggerPathToPass    = parameters.getParameter<string>("muonTriggerPathToPass");
00056   theTriggerResultsCollection = parameters.getParameter<InputTag>("triggerResultsCollection");
00057   theMuonCollectionLabel      = parameters.getParameter<InputTag>("muonCollection");
00058   theElectronCollectionLabel  = parameters.getParameter<InputTag>("electronCollection");
00059   theCaloJetCollectionLabel   = parameters.getParameter<InputTag>("caloJetCollection");
00060   theCaloMETCollectionLabel   = parameters.getParameter<InputTag>("caloMETCollection");
00061   
00062   // just to initialize
00063   isValidHltConfig_ = false;
00064 
00065 }
00066 
00067 EwkDQM::~EwkDQM() { 
00068 }
00069 
00070 
00071 void EwkDQM::beginJob() {
00072 
00073   logTraceName = "EwkAnalyzer";
00074 
00075   LogTrace(logTraceName)<<"Parameters initialization";
00076   theDbe = Service<DQMStore>().operator->();
00077   theDbe->setCurrentFolder("Physics/EwkDQM");  // Use folder with name of PAG
00078 
00079   const float pi = 3.14159265;
00080 
00081   // Keep the number of plots and number of bins to a minimum!
00082   h_vertex_number = theDbe->book1D("h_vertex_number", "Number of event vertices in collection", 10,-0.5,   9.5 );
00083   h_vertex_chi2  = theDbe->book1D("h_vertex_chi2" , "Event Vertex #chi^{2}/n.d.o.f."          , 20, 0.0,   2.0 );
00084   h_vertex_numTrks = theDbe->book1D("h_vertex_numTrks", "Event Vertex, number of tracks"     , 20, -0.5,  59.5 );
00085   h_vertex_sumTrks = theDbe->book1D("h_vertex_sumTrks", "Event Vertex, sum of track pt"      , 20,  0.0, 100.0 );
00086   h_vertex_d0    = theDbe->book1D("h_vertex_d0"   , "Event Vertex d0"                        , 20,  0.0,   0.05);
00087   h_mumu_invMass = theDbe->book1D("h_mumu_invMass", "#mu#mu Invariant Mass;InvMass (GeV)"    , 20, 40.0, 140.0 );
00088   h_ee_invMass   = theDbe->book1D("h_ee_invMass",   "ee Invariant Mass;InvMass (Gev)"        , 20, 40.0, 140.0 );
00089   h_jet_et       = theDbe->book1D("h_jet_et",       "Jet with highest E_{T} (from "+theCaloJetCollectionLabel.label()+");E_{T}(1^{st} jet) (GeV)",    20, 0., 200.0);
00090   h_jet2_et      = theDbe->book1D("h_jet2_et",      "Jet with 2^{nd} highest E_{T} (from "+theCaloJetCollectionLabel.label()+");E_{T}(2^{nd} jet) (GeV)",    20, 0., 200.0);
00091   h_jet_count    = theDbe->book1D("h_jet_count",    "Number of "+theCaloJetCollectionLabel.label()+" (E_{T} > 15 GeV);Number of Jets", 8, -0.5, 7.5);
00092   h_e1_et        = theDbe->book1D("h_e1_et",  "E_{T} of Leading Electron;E_{T} (GeV)"        , 20,  0.0 , 100.0);
00093   h_e2_et        = theDbe->book1D("h_e2_et",  "E_{T} of Second Electron;E_{T} (GeV)"         , 20,  0.0 , 100.0);
00094   h_e1_eta       = theDbe->book1D("h_e1_eta", "#eta of Leading Electron;#eta"                , 20, -4.0 , 4.0);
00095   h_e2_eta       = theDbe->book1D("h_e2_eta", "#eta of Second Electron;#eta"                 , 20, -4.0 , 4.0);
00096   h_e1_phi       = theDbe->book1D("h_e1_phi", "#phi of Leading Electron;#phi"                , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
00097   h_e2_phi       = theDbe->book1D("h_e2_phi", "#phi of Second Electron;#phi"                 , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
00098   h_m1_pt        = theDbe->book1D("h_m1_pt",  "p_{T} of Leading Muon;p_{T}(1^{st} #mu) (GeV)", 20,  0.0 , 100.0);
00099   h_m2_pt        = theDbe->book1D("h_m2_pt",  "p_{T} of Second Muon;p_{T}(2^{nd} #mu) (GeV)" , 20,  0.0 , 100.0);
00100   h_m1_eta       = theDbe->book1D("h_m1_eta", "#eta of Leading Muon;#eta(1^{st} #mu)"        , 20, -4.0 , 4.0);
00101   h_m2_eta       = theDbe->book1D("h_m2_eta", "#eta of Second Muon;#eta(2^{nd} #mu)"         , 20, -4.0 , 4.0);
00102   h_m1_phi       = theDbe->book1D("h_m1_phi", "#phi of Leading Muon;#phi(1^{st} #mu)"        , 20, (-1.-1./10.)*pi, (1.+1./10.)*pi);
00103   h_m2_phi       = theDbe->book1D("h_m2_phi", "#phi of Second Muon;#phi(2^{nd} #mu)"         , 20, (-1.-1./10.)*pi, (1.+1./10.)*pi);
00104 //  h_t1_et          = theDbe->book1D("h_t1_et",           "E_{T} of Leading Tau;E_{T} (GeV)" , 20, 0.0 , 100.0);
00105 //  h_t1_eta         = theDbe->book1D("h_t1_eta",          "#eta of Leading Tau;#eta"               , 20, -4.0, 4.0);
00106 //  h_t1_phi         = theDbe->book1D("h_t1_phi",          "#phi of Leading Tau;#phi"               , 20, -4.0, 4.0);
00107   h_met          = theDbe->book1D("h_met",        "Missing E_{T}; GeV"                       , 20,  0.0 , 100);
00108   h_met_phi      = theDbe->book1D("h_met_phi",    "Missing E_{T} #phi;#phi(MET)"             , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
00109   h_e_invWMass   = theDbe->book1D("h_e_invWMass", "W-> e #nu Transverse Mass;M_{T} (GeV)"    , 20,  0.0, 140.0); 
00110   h_m_invWMass   = theDbe->book1D("h_m_invWMass", "W-> #mu #nu Transverse Mass;M_{T} (GeV)"  , 20,  0.0, 140.0); 
00111 }
00112 
00113 
00114 
00118 void EwkDQM::beginRun( const edm::Run& theRun, const edm::EventSetup& theSetup ) {
00119   
00120   // passed as parameter to HLTConfigProvider::init(), not yet used
00121   bool isConfigChanged = false;
00122   
00123   // isValidHltConfig_ used to short-circuit analyze() in case of problems
00124   const std::string hltProcessName( theTriggerResultsCollection.process() );
00125   isValidHltConfig_ = hltConfigProvider_.init( theRun, theSetup, hltProcessName, isConfigChanged );
00126 
00127 }
00128 
00129 
00130 void EwkDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
00131 
00132   // short-circuit if hlt problems
00133   if( ! isValidHltConfig_ ) return;
00134 
00135 
00136   LogTrace(logTraceName)<<"Analysis of event # ";
00137   // Did it pass certain HLT path?
00138   Handle<TriggerResults> HLTresults;
00139   iEvent.getByLabel(theTriggerResultsCollection, HLTresults); 
00140   if ( !HLTresults.isValid() ) return;
00141 
00142   //unsigned int triggerIndex_elec = hltConfig.triggerIndex(theElecTriggerPathToPass);
00143   //unsigned int triggerIndex_muon = hltConfig.triggerIndex(theMuonTriggerPathToPass);
00144   bool passed_electron_HLT = true;
00145   bool passed_muon_HLT     = true;
00146   //if (triggerIndex_elec < HLTresults->size()) passed_electron_HLT = HLTresults->accept(triggerIndex_elec);
00147   //if (triggerIndex_muon < HLTresults->size()) passed_muon_HLT     = HLTresults->accept(triggerIndex_muon);
00148   //if ( !(passed_electron_HLT || passed_muon_HLT) ) return;
00149 
00151   //Vertex information
00152   Handle<VertexCollection> vertexHandle;
00153   iEvent.getByLabel("offlinePrimaryVertices", vertexHandle);
00154   if ( !vertexHandle.isValid() ) return;
00155   VertexCollection vertexCollection = *(vertexHandle.product());
00156   int vertex_number     = vertexCollection.size();
00157   VertexCollection::const_iterator v = vertexCollection.begin();
00158   double vertex_chi2    = v->normalizedChi2(); //v->chi2();
00159   double vertex_d0      = sqrt(v->x()*v->x()+v->y()*v->y());
00160   //double vertex_ndof    = v->ndof();cout << "ndof="<<vertex_ndof<<endl;
00161   double vertex_numTrks = v->tracksSize();
00162   double vertex_sumTrks = 0.0;
00163   for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack!=v->tracks_end(); vertex_curTrack++) {
00164     vertex_sumTrks += (*vertex_curTrack)->pt();
00165   }
00166 
00168   //Missing ET
00169   Handle<CaloMETCollection> caloMETCollection;
00170   iEvent.getByLabel(theCaloMETCollectionLabel, caloMETCollection);
00171   if ( !caloMETCollection.isValid() ) return;
00172   float missing_et = caloMETCollection->begin()->et();
00173   float met_phi    = caloMETCollection->begin()->phi();
00174 
00175 
00177   // grab "gaussian sum fitting" electrons
00178   Handle<GsfElectronCollection> electronCollection;
00179   iEvent.getByLabel(theElectronCollectionLabel, electronCollection);
00180   if ( !electronCollection.isValid() ) return;
00181 
00182   // Find the highest and 2nd highest electron
00183   float electron_et   = -8.0;
00184   float electron_eta  = -8.0;
00185   float electron_phi  = -8.0;
00186   float electron2_et  = -9.0;
00187   float electron2_eta = -9.0;
00188   float electron2_phi = -9.0;
00189   float ee_invMass    = -9.0;
00190   TLorentzVector e1, e2;
00191 
00192   // If it passed electron HLT and the collection was found, find electrons near Z mass
00193   if( passed_electron_HLT ) {
00194 
00195     for (reco::GsfElectronCollection::const_iterator recoElectron=electronCollection->begin(); recoElectron!=electronCollection->end(); recoElectron++){
00196 
00197       // Require electron to pass some basic cuts
00198       if ( recoElectron->et() < 20 || fabs(recoElectron->eta())>2.5 ) continue;
00199 
00200       // Tighter electron cuts
00201       if ( recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 || 
00202            recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 || 
00203            recoElectron->sigmaIetaIeta() > 0.027 ) continue;
00204 
00205       if (recoElectron->et() > electron_et){
00206         electron2_et  = electron_et;  // 2nd highest gets values from current highest
00207         electron2_eta = electron_eta;
00208         electron2_phi = electron_phi;
00209         electron_et   = recoElectron->et();  // 1st highest gets values from new highest
00210         electron_eta  = recoElectron->eta();
00211         electron_phi  = recoElectron->phi();
00212         e1 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
00213       } else if (recoElectron->et() > electron2_et) {
00214         electron2_et  = recoElectron->et();
00215         electron2_eta = recoElectron->eta();
00216         electron2_phi = recoElectron->phi();
00217         e2 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
00218       }
00219     } // end of loop over electrons
00220     if (electron2_et>0.0) {
00221       TLorentzVector pair=e1+e2;
00222       ee_invMass = pair.M();
00223     }
00224   } // end of "are electrons valid"
00226 
00227 
00228 
00230   // Take the STA muon container
00231   Handle<MuonCollection> muonCollection;
00232   iEvent.getByLabel(theMuonCollectionLabel,muonCollection);
00233   if ( !muonCollection.isValid() ) return;
00234 
00235   // Find the highest pt muons
00236   float mm_invMass = -9.0;
00237   float muon_pt   = -9.0;
00238   float muon_eta  = -9.0;
00239   float muon_phi  = -9.0;
00240   float muon2_pt  = -9.0;
00241   float muon2_eta = -9.0;
00242   float muon2_phi = -9.0;
00243   TLorentzVector m1, m2;
00244 
00245   if( passed_muon_HLT ) {
00246     for (reco::MuonCollection::const_iterator recoMuon=muonCollection->begin(); recoMuon!=muonCollection->end(); recoMuon++){
00247 
00248       // Require muon to pass some basic cuts
00249       if ( recoMuon->pt() < 20 || !recoMuon->isGlobalMuon() ) continue;
00250       // Some tighter muon cuts
00251       if ( recoMuon->globalTrack()->normalizedChi2() > 10 ) continue;
00252 
00253       if (recoMuon->pt() > muon_pt){
00254         muon2_pt  = muon_pt;  // 2nd highest gets values from current highest    
00255         muon2_eta = muon_eta;
00256         muon2_phi = muon_phi;
00257         muon_pt   = recoMuon->pt();  // 1st highest gets values from new highest
00258         muon_eta  = recoMuon->eta();
00259         muon_phi  = recoMuon->phi();
00260         m1 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
00261       } else if (recoMuon->pt() > muon2_pt) {
00262         muon2_pt  = recoMuon->pt();
00263         muon2_eta = recoMuon->eta();
00264         muon2_phi = recoMuon->phi();
00265         m2 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
00266       }
00267     }
00268   }
00269   if (muon2_pt>0.0) {
00270     TLorentzVector pair=m1+m2;
00271     mm_invMass = pair.M();
00272   }
00274 
00275 
00277   // Find the highest et jet
00278   Handle<CaloJetCollection> caloJetCollection;
00279   iEvent.getByLabel (theCaloJetCollectionLabel,caloJetCollection);
00280   if ( !caloJetCollection.isValid() ) return;
00281 
00282   float jet_et    = -8.0;
00283   float jet_eta   = -8.0;
00284   float jet_phi   = -8.0;
00285   int   jet_count = 0;
00286   float jet2_et   = -9.0;
00287   float jet2_eta  = -9.0;
00288   float jet2_phi  = -9.0;
00289   for (CaloJetCollection::const_iterator i_calojet = caloJetCollection->begin(); i_calojet != caloJetCollection->end(); i_calojet++) {
00290 
00291     float jet_current_et = i_calojet->et();
00292 
00293     // if it overlaps with electron, it is not a jet
00294     if ( electron_et>0.0 && fabs(i_calojet->eta()-electron_eta ) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron_phi ) < 0.2) continue;
00295     if ( electron2_et>0.0&& fabs(i_calojet->eta()-electron2_eta) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron2_phi) < 0.2) continue;
00296 
00297     // if it has too low Et, throw away
00298     if (jet_current_et < 15) continue;
00299 
00300     jet_count++;
00301     if (jet_current_et > jet_et) {
00302       jet2_et  = jet_et;  // 2nd highest jet get's et from current highest
00303       jet2_eta = jet_eta;
00304       jet2_phi = jet_phi;
00305       jet_et   = i_calojet->et(); // current highest jet gets et from the new highest
00306       jet_eta  = i_calojet->eta();
00307       jet_phi  = i_calojet->phi();
00308     } else if (jet_current_et > jet2_et) {
00309       jet2_et  = i_calojet->et();
00310       jet2_eta = i_calojet->eta();
00311       jet2_phi = i_calojet->phi();
00312     }
00313   }
00315 
00316 
00317 
00319   //                 Fill Histograms                                            //
00321 
00322   bool fill_e1  = false;
00323   bool fill_e2  = false;
00324   bool fill_m1  = false;
00325   bool fill_m2  = false;
00326   bool fill_met = false;
00327 
00328   // Was Z->ee found?
00329   if (ee_invMass>0.0) {
00330     h_ee_invMass ->Fill(ee_invMass);
00331     fill_e1 = true;
00332     fill_e2 = true;
00333   }
00334 
00335   // Was Z->mu mu found?
00336   if (mm_invMass > 0.0) {
00337     h_mumu_invMass->Fill(mm_invMass);
00338     fill_m1 = true;
00339     fill_m2 = true;
00340   }
00341 
00342   // Was W->e nu found?
00343   if (electron_et>0.0&&missing_et>20.0) {
00344     float dphiW = fabs(met_phi-electron_phi);
00345     float W_mt_e = sqrt(2*missing_et*electron_et*(1-cos(dphiW)));
00346     h_e_invWMass ->Fill(W_mt_e);
00347     fill_e1  = true;
00348     fill_met = true;
00349   }
00350 
00351   // Was W->mu nu found?
00352   if (muon_pt>0.0&&missing_et>20.0) {
00353     float dphiW = fabs(met_phi-muon_phi);
00354     float W_mt_m = sqrt(2*missing_et*muon_pt*(1-cos(dphiW)));
00355     h_m_invWMass ->Fill(W_mt_m);
00356     fill_m1  = true;
00357     fill_met = true;
00358   }
00359 
00360   if (jet_et>0.0) {
00361     h_jet_et   ->Fill(jet_et);
00362     h_jet_count->Fill(jet_count);
00363   }
00364 
00365   if (fill_e1 || fill_m1) {
00366     h_vertex_number->Fill(vertex_number);
00367     h_vertex_chi2->Fill(vertex_chi2);
00368     h_vertex_d0  ->Fill(vertex_d0);
00369     h_vertex_numTrks->Fill(vertex_numTrks);
00370     h_vertex_sumTrks->Fill(vertex_sumTrks);
00371   }
00372 
00373   if (fill_e1) {
00374     h_e1_et      ->Fill(electron_et);
00375     h_e1_eta     ->Fill(electron_eta);
00376     h_e1_phi     ->Fill(electron_phi);
00377   }
00378   if (fill_e2) {
00379     h_e2_et      ->Fill(electron2_et);
00380     h_e2_eta     ->Fill(electron2_eta);
00381     h_e2_phi     ->Fill(electron2_phi);
00382   }
00383   if (fill_m1) {
00384     h_m1_pt      ->Fill(muon_pt);
00385     h_m1_eta     ->Fill(muon_eta);
00386     h_m1_phi     ->Fill(muon_phi);
00387   }
00388   if (fill_m2) {
00389     h_m2_pt      ->Fill(muon2_pt);
00390     h_m2_eta     ->Fill(muon2_eta);
00391     h_m2_phi     ->Fill(muon2_phi);
00392   }
00393   if (fill_met) {
00394     h_met        ->Fill(missing_et);
00395     h_met_phi    ->Fill(met_phi);
00396   }
00398 }
00399 
00400 
00401 void EwkDQM::endJob(void) {}
00402 
00403 
00404 // This always returns only a positive deltaPhi
00405 double EwkDQM::calcDeltaPhi(double phi1, double phi2) {
00406 
00407   double deltaPhi = phi1 - phi2;
00408 
00409   if (deltaPhi < 0) deltaPhi = -deltaPhi;
00410 
00411   if (deltaPhi > 3.1415926) {
00412     deltaPhi = 2 * 3.1415926 - deltaPhi;
00413   }
00414 
00415   return deltaPhi;
00416 }
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/DQM/Physics/src/EwkDQM.cc
