CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EwkDQM.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * $Date: 2010/03/12 10:25:45 $
5  * $Revision: 1.15 $
6  * \author Michael B. Anderson, University of Wisconsin-Madison
7  * \author Will Parker, University of Wisconsin-Madison
8  */
9 
10 #include "DQM/Physics/src/EwkDQM.h"
15 
19 
21 
23 
24 // Physics Objects
37 
38 
40 #include "TLorentzVector.h"
41 
42 #include <vector>
43 
44 #include <string>
45 #include <cmath>
46 using namespace std;
47 using namespace edm;
48 using namespace reco;
49 
50 
51 
53  // Get parameters from configuration file
54  theElecTriggerPathToPass = parameters.getParameter<string>("elecTriggerPathToPass");
55  theMuonTriggerPathToPass = parameters.getParameter<string>("muonTriggerPathToPass");
56  theTriggerResultsCollection = parameters.getParameter<InputTag>("triggerResultsCollection");
57  theMuonCollectionLabel = parameters.getParameter<InputTag>("muonCollection");
58  theElectronCollectionLabel = parameters.getParameter<InputTag>("electronCollection");
59  theCaloJetCollectionLabel = parameters.getParameter<InputTag>("caloJetCollection");
60  theCaloMETCollectionLabel = parameters.getParameter<InputTag>("caloMETCollection");
61 
62  // just to initialize
63  isValidHltConfig_ = false;
64 
65 }
66 
68 }
69 
70 
72 
73  logTraceName = "EwkAnalyzer";
74 
75  LogTrace(logTraceName)<<"Parameters initialization";
76  theDbe = Service<DQMStore>().operator->();
77  theDbe->setCurrentFolder("Physics/EwkDQM"); // Use folder with name of PAG
78 
79  const float pi = 3.14159265;
80 
81  // Keep the number of plots and number of bins to a minimum!
82  h_vertex_number = theDbe->book1D("h_vertex_number", "Number of event vertices in collection", 10,-0.5, 9.5 );
83  h_vertex_chi2 = theDbe->book1D("h_vertex_chi2" , "Event Vertex #chi^{2}/n.d.o.f." , 20, 0.0, 2.0 );
84  h_vertex_numTrks = theDbe->book1D("h_vertex_numTrks", "Event Vertex, number of tracks" , 20, -0.5, 59.5 );
85  h_vertex_sumTrks = theDbe->book1D("h_vertex_sumTrks", "Event Vertex, sum of track pt" , 20, 0.0, 100.0 );
86  h_vertex_d0 = theDbe->book1D("h_vertex_d0" , "Event Vertex d0" , 20, 0.0, 0.05);
87  h_mumu_invMass = theDbe->book1D("h_mumu_invMass", "#mu#mu Invariant Mass;InvMass (GeV)" , 20, 40.0, 140.0 );
88  h_ee_invMass = theDbe->book1D("h_ee_invMass", "ee Invariant Mass;InvMass (Gev)" , 20, 40.0, 140.0 );
89  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);
90  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);
91  h_jet_count = theDbe->book1D("h_jet_count", "Number of "+theCaloJetCollectionLabel.label()+" (E_{T} > 15 GeV);Number of Jets", 8, -0.5, 7.5);
92  h_e1_et = theDbe->book1D("h_e1_et", "E_{T} of Leading Electron;E_{T} (GeV)" , 20, 0.0 , 100.0);
93  h_e2_et = theDbe->book1D("h_e2_et", "E_{T} of Second Electron;E_{T} (GeV)" , 20, 0.0 , 100.0);
94  h_e1_eta = theDbe->book1D("h_e1_eta", "#eta of Leading Electron;#eta" , 20, -4.0 , 4.0);
95  h_e2_eta = theDbe->book1D("h_e2_eta", "#eta of Second Electron;#eta" , 20, -4.0 , 4.0);
96  h_e1_phi = theDbe->book1D("h_e1_phi", "#phi of Leading Electron;#phi" , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
97  h_e2_phi = theDbe->book1D("h_e2_phi", "#phi of Second Electron;#phi" , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
98  h_m1_pt = theDbe->book1D("h_m1_pt", "p_{T} of Leading Muon;p_{T}(1^{st} #mu) (GeV)", 20, 0.0 , 100.0);
99  h_m2_pt = theDbe->book1D("h_m2_pt", "p_{T} of Second Muon;p_{T}(2^{nd} #mu) (GeV)" , 20, 0.0 , 100.0);
100  h_m1_eta = theDbe->book1D("h_m1_eta", "#eta of Leading Muon;#eta(1^{st} #mu)" , 20, -4.0 , 4.0);
101  h_m2_eta = theDbe->book1D("h_m2_eta", "#eta of Second Muon;#eta(2^{nd} #mu)" , 20, -4.0 , 4.0);
102  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);
103  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);
104 // h_t1_et = theDbe->book1D("h_t1_et", "E_{T} of Leading Tau;E_{T} (GeV)" , 20, 0.0 , 100.0);
105 // h_t1_eta = theDbe->book1D("h_t1_eta", "#eta of Leading Tau;#eta" , 20, -4.0, 4.0);
106 // h_t1_phi = theDbe->book1D("h_t1_phi", "#phi of Leading Tau;#phi" , 20, -4.0, 4.0);
107  h_met = theDbe->book1D("h_met", "Missing E_{T}; GeV" , 20, 0.0 , 100);
108  h_met_phi = theDbe->book1D("h_met_phi", "Missing E_{T} #phi;#phi(MET)" , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
109  h_e_invWMass = theDbe->book1D("h_e_invWMass", "W-> e #nu Transverse Mass;M_{T} (GeV)" , 20, 0.0, 140.0);
110  h_m_invWMass = theDbe->book1D("h_m_invWMass", "W-> #mu #nu Transverse Mass;M_{T} (GeV)" , 20, 0.0, 140.0);
111 }
112 
113 
114 
118 void EwkDQM::beginRun( const edm::Run& theRun, const edm::EventSetup& theSetup ) {
119 
120  // passed as parameter to HLTConfigProvider::init(), not yet used
121  bool isConfigChanged = false;
122 
123  // isValidHltConfig_ used to short-circuit analyze() in case of problems
124  const std::string hltProcessName( theTriggerResultsCollection.process() );
125  isValidHltConfig_ = hltConfigProvider_.init( theRun, theSetup, hltProcessName, isConfigChanged );
126 
127 }
128 
129 
130 void EwkDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
131 
132  // short-circuit if hlt problems
133  if( ! isValidHltConfig_ ) return;
134 
135 
136  LogTrace(logTraceName)<<"Analysis of event # ";
137  // Did it pass certain HLT path?
138  Handle<TriggerResults> HLTresults;
139  iEvent.getByLabel(theTriggerResultsCollection, HLTresults);
140  if ( !HLTresults.isValid() ) return;
141 
142  //unsigned int triggerIndex_elec = hltConfig.triggerIndex(theElecTriggerPathToPass);
143  //unsigned int triggerIndex_muon = hltConfig.triggerIndex(theMuonTriggerPathToPass);
144  bool passed_electron_HLT = true;
145  bool passed_muon_HLT = true;
146  //if (triggerIndex_elec < HLTresults->size()) passed_electron_HLT = HLTresults->accept(triggerIndex_elec);
147  //if (triggerIndex_muon < HLTresults->size()) passed_muon_HLT = HLTresults->accept(triggerIndex_muon);
148  //if ( !(passed_electron_HLT || passed_muon_HLT) ) return;
149 
151  //Vertex information
152  Handle<VertexCollection> vertexHandle;
153  iEvent.getByLabel("offlinePrimaryVertices", vertexHandle);
154  if ( !vertexHandle.isValid() ) return;
155  VertexCollection vertexCollection = *(vertexHandle.product());
156  int vertex_number = vertexCollection.size();
157  VertexCollection::const_iterator v = vertexCollection.begin();
158  double vertex_chi2 = v->normalizedChi2(); //v->chi2();
159  double vertex_d0 = sqrt(v->x()*v->x()+v->y()*v->y());
160  //double vertex_ndof = v->ndof();cout << "ndof="<<vertex_ndof<<endl;
161  double vertex_numTrks = v->tracksSize();
162  double vertex_sumTrks = 0.0;
163  for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack!=v->tracks_end(); vertex_curTrack++) {
164  vertex_sumTrks += (*vertex_curTrack)->pt();
165  }
166 
168  //Missing ET
169  Handle<CaloMETCollection> caloMETCollection;
170  iEvent.getByLabel(theCaloMETCollectionLabel, caloMETCollection);
171  if ( !caloMETCollection.isValid() ) return;
172  float missing_et = caloMETCollection->begin()->et();
173  float met_phi = caloMETCollection->begin()->phi();
174 
175 
177  // grab "gaussian sum fitting" electrons
178  Handle<GsfElectronCollection> electronCollection;
179  iEvent.getByLabel(theElectronCollectionLabel, electronCollection);
180  if ( !electronCollection.isValid() ) return;
181 
182  // Find the highest and 2nd highest electron
183  float electron_et = -8.0;
184  float electron_eta = -8.0;
185  float electron_phi = -8.0;
186  float electron2_et = -9.0;
187  float electron2_eta = -9.0;
188  float electron2_phi = -9.0;
189  float ee_invMass = -9.0;
190  TLorentzVector e1, e2;
191 
192  // If it passed electron HLT and the collection was found, find electrons near Z mass
193  if( passed_electron_HLT ) {
194 
195  for (reco::GsfElectronCollection::const_iterator recoElectron=electronCollection->begin(); recoElectron!=electronCollection->end(); recoElectron++){
196 
197  // Require electron to pass some basic cuts
198  if ( recoElectron->et() < 20 || fabs(recoElectron->eta())>2.5 ) continue;
199 
200  // Tighter electron cuts
201  if ( recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 ||
202  recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 ||
203  recoElectron->sigmaIetaIeta() > 0.027 ) continue;
204 
205  if (recoElectron->et() > electron_et){
206  electron2_et = electron_et; // 2nd highest gets values from current highest
207  electron2_eta = electron_eta;
208  electron2_phi = electron_phi;
209  electron_et = recoElectron->et(); // 1st highest gets values from new highest
210  electron_eta = recoElectron->eta();
211  electron_phi = recoElectron->phi();
212  e1 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
213  } else if (recoElectron->et() > electron2_et) {
214  electron2_et = recoElectron->et();
215  electron2_eta = recoElectron->eta();
216  electron2_phi = recoElectron->phi();
217  e2 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
218  }
219  } // end of loop over electrons
220  if (electron2_et>0.0) {
221  TLorentzVector pair=e1+e2;
222  ee_invMass = pair.M();
223  }
224  } // end of "are electrons valid"
226 
227 
228 
230  // Take the STA muon container
231  Handle<MuonCollection> muonCollection;
232  iEvent.getByLabel(theMuonCollectionLabel,muonCollection);
233  if ( !muonCollection.isValid() ) return;
234 
235  // Find the highest pt muons
236  float mm_invMass = -9.0;
237  float muon_pt = -9.0;
238  float muon_eta = -9.0;
239  float muon_phi = -9.0;
240  float muon2_pt = -9.0;
241  float muon2_eta = -9.0;
242  float muon2_phi = -9.0;
243  TLorentzVector m1, m2;
244 
245  if( passed_muon_HLT ) {
246  for (reco::MuonCollection::const_iterator recoMuon=muonCollection->begin(); recoMuon!=muonCollection->end(); recoMuon++){
247 
248  // Require muon to pass some basic cuts
249  if ( recoMuon->pt() < 20 || !recoMuon->isGlobalMuon() ) continue;
250  // Some tighter muon cuts
251  if ( recoMuon->globalTrack()->normalizedChi2() > 10 ) continue;
252 
253  if (recoMuon->pt() > muon_pt){
254  muon2_pt = muon_pt; // 2nd highest gets values from current highest
255  muon2_eta = muon_eta;
256  muon2_phi = muon_phi;
257  muon_pt = recoMuon->pt(); // 1st highest gets values from new highest
258  muon_eta = recoMuon->eta();
259  muon_phi = recoMuon->phi();
260  m1 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
261  } else if (recoMuon->pt() > muon2_pt) {
262  muon2_pt = recoMuon->pt();
263  muon2_eta = recoMuon->eta();
264  muon2_phi = recoMuon->phi();
265  m2 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
266  }
267  }
268  }
269  if (muon2_pt>0.0) {
270  TLorentzVector pair=m1+m2;
271  mm_invMass = pair.M();
272  }
274 
275 
277  // Find the highest et jet
278  Handle<CaloJetCollection> caloJetCollection;
279  iEvent.getByLabel (theCaloJetCollectionLabel,caloJetCollection);
280  if ( !caloJetCollection.isValid() ) return;
281 
282  float jet_et = -8.0;
283  float jet_eta = -8.0;
284  float jet_phi = -8.0;
285  int jet_count = 0;
286  float jet2_et = -9.0;
287  float jet2_eta = -9.0;
288  float jet2_phi = -9.0;
289  for (CaloJetCollection::const_iterator i_calojet = caloJetCollection->begin(); i_calojet != caloJetCollection->end(); i_calojet++) {
290 
291  float jet_current_et = i_calojet->et();
292 
293  // if it overlaps with electron, it is not a jet
294  if ( electron_et>0.0 && fabs(i_calojet->eta()-electron_eta ) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron_phi ) < 0.2) continue;
295  if ( electron2_et>0.0&& fabs(i_calojet->eta()-electron2_eta) < 0.2 && calcDeltaPhi(i_calojet->phi(), electron2_phi) < 0.2) continue;
296 
297  // if it has too low Et, throw away
298  if (jet_current_et < 15) continue;
299 
300  jet_count++;
301  if (jet_current_et > jet_et) {
302  jet2_et = jet_et; // 2nd highest jet get's et from current highest
303  jet2_eta = jet_eta;
304  jet2_phi = jet_phi;
305  jet_et = i_calojet->et(); // current highest jet gets et from the new highest
306  jet_eta = i_calojet->eta();
307  jet_phi = i_calojet->phi();
308  } else if (jet_current_et > jet2_et) {
309  jet2_et = i_calojet->et();
310  jet2_eta = i_calojet->eta();
311  jet2_phi = i_calojet->phi();
312  }
313  }
315 
316 
317 
319  // Fill Histograms //
321 
322  bool fill_e1 = false;
323  bool fill_e2 = false;
324  bool fill_m1 = false;
325  bool fill_m2 = false;
326  bool fill_met = false;
327 
328  // Was Z->ee found?
329  if (ee_invMass>0.0) {
330  h_ee_invMass ->Fill(ee_invMass);
331  fill_e1 = true;
332  fill_e2 = true;
333  }
334 
335  // Was Z->mu mu found?
336  if (mm_invMass > 0.0) {
337  h_mumu_invMass->Fill(mm_invMass);
338  fill_m1 = true;
339  fill_m2 = true;
340  }
341 
342  // Was W->e nu found?
343  if (electron_et>0.0&&missing_et>20.0) {
344  float dphiW = fabs(met_phi-electron_phi);
345  float W_mt_e = sqrt(2*missing_et*electron_et*(1-cos(dphiW)));
346  h_e_invWMass ->Fill(W_mt_e);
347  fill_e1 = true;
348  fill_met = true;
349  }
350 
351  // Was W->mu nu found?
352  if (muon_pt>0.0&&missing_et>20.0) {
353  float dphiW = fabs(met_phi-muon_phi);
354  float W_mt_m = sqrt(2*missing_et*muon_pt*(1-cos(dphiW)));
355  h_m_invWMass ->Fill(W_mt_m);
356  fill_m1 = true;
357  fill_met = true;
358  }
359 
360  if (jet_et>0.0) {
361  h_jet_et ->Fill(jet_et);
362  h_jet_count->Fill(jet_count);
363  }
364 
365  if (fill_e1 || fill_m1) {
366  h_vertex_number->Fill(vertex_number);
367  h_vertex_chi2->Fill(vertex_chi2);
368  h_vertex_d0 ->Fill(vertex_d0);
369  h_vertex_numTrks->Fill(vertex_numTrks);
370  h_vertex_sumTrks->Fill(vertex_sumTrks);
371  }
372 
373  if (fill_e1) {
374  h_e1_et ->Fill(electron_et);
375  h_e1_eta ->Fill(electron_eta);
376  h_e1_phi ->Fill(electron_phi);
377  }
378  if (fill_e2) {
379  h_e2_et ->Fill(electron2_et);
380  h_e2_eta ->Fill(electron2_eta);
381  h_e2_phi ->Fill(electron2_phi);
382  }
383  if (fill_m1) {
384  h_m1_pt ->Fill(muon_pt);
385  h_m1_eta ->Fill(muon_eta);
386  h_m1_phi ->Fill(muon_phi);
387  }
388  if (fill_m2) {
389  h_m2_pt ->Fill(muon2_pt);
390  h_m2_eta ->Fill(muon2_eta);
391  h_m2_phi ->Fill(muon2_phi);
392  }
393  if (fill_met) {
394  h_met ->Fill(missing_et);
395  h_met_phi ->Fill(met_phi);
396  }
398 }
399 
400 
401 void EwkDQM::endJob(void) {}
402 
403 
404 // This always returns only a positive deltaPhi
405 double EwkDQM::calcDeltaPhi(double phi1, double phi2) {
406 
407  double deltaPhi = phi1 - phi2;
408 
409  if (deltaPhi < 0) deltaPhi = -deltaPhi;
410 
411  if (deltaPhi > 3.1415926) {
412  deltaPhi = 2 * 3.1415926 - deltaPhi;
413  }
414 
415  return deltaPhi;
416 }
void beginRun(const edm::Run &, const edm::EventSetup &)
Definition: EwkDQM.cc:118
T getParameter(std::string const &) const
double calcDeltaPhi(double phi1, double phi2)
Definition: EwkTauDQM.cc:972
dictionary parameters
Definition: Parameters.py:2
EwkDQM(const edm::ParameterSet &)
Constructor.
Definition: EwkDQM.cc:52
double deltaPhi(float phi1, float phi2)
Definition: VectorUtil.h:30
void endJob(void)
Save the histos.
Definition: EwkDQM.cc:401
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
void analyze(const edm::Event &, const edm::EventSetup &)
Get the analysis.
Definition: EwkDQM.cc:130
tuple vertexCollection
int iEvent
Definition: GenABIO.cc:243
T sqrt(T t)
Definition: SSEVec.h:28
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
void beginJob()
Inizialize parameters for histo binning.
Definition: EwkDQM.cc:71
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
#define LogTrace(id)
virtual ~EwkDQM()
Destructor.
Definition: EwkDQM.cc:67
T const * product() const
Definition: Handle.h:74
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:38
double calcDeltaPhi(double phi1, double phi2)
Definition: EwkDQM.cc:405
double pi
mathSSE::Vec4< T > v
Definition: Run.h:32