CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MCEmbeddingValidationAnalyzer.cc
Go to the documentation of this file.
2 
5 
21 
23 
25 
26 #include <Math/VectorUtil.h>
27 #include <TMath.h>
28 
29 #include <sstream>
30 
31 // Copied from TauAnalysis/CandidateTools/interface/svFitAuxFunctions.cc. Note that we cannot use
32 // this package directly, since it is not part of CMSSW, but we are.
33 namespace SVfit_namespace
34 {
35  inline double square(double x)
36  {
37  return x*x;
38  }
39 
41  {
42  double p_x = p.x();
43  double p_y = p.y();
44  double p_z = p.z();
45  double mag2 = square(p_x) + square(p_y) + square(p_z);
46  if ( mag2 <= 0. ) return p;
47  double mag = TMath::Sqrt(mag2);
48  return reco::Candidate::Vector(p_x/mag, p_y/mag, p_z/mag);
49  }
50 
52  {
53  return (p1.x()*p2.x() + p1.y()*p2.y() + p1.z()*p2.z());
54  }
55 
57  {
58  double p3_x = p1.y()*p2.z() - p1.z()*p2.y();
59  double p3_y = p1.z()*p2.x() - p1.x()*p2.z();
60  double p3_z = p1.x()*p2.y() - p1.y()*p2.x();
61  return reco::Candidate::Vector(p3_x, p3_y, p3_z);
62  }
63 
65  {
66  reco::Candidate::Vector u_z = normalize(reco::Candidate::Vector(visP4.px(), visP4.py(), visP4.pz()));
69 
70  reco::Candidate::Vector p3Mother_unit = normalize(reco::Candidate::Vector(motherP4.px(), motherP4.py(), motherP4.pz()));
71 
72  double phi_lab = TMath::ATan2(compScalarProduct(p3Mother_unit, u_y), compScalarProduct(p3Mother_unit, u_x));
73  return phi_lab;
74  }
75 }
76 
77 /*EGammaMvaEleEstimator* MCEmbeddingValidationAnalyzer::electronDistributionExtra::fMVA_ = 0;
78 bool MCEmbeddingValidationAnalyzer::electronDistributionExtra::fMVA_isInitialized_ = false;*/
79 
81  : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
82  srcReplacedMuons_(cfg.getParameter<edm::InputTag>("srcReplacedMuons")),
83  srcRecMuons_(cfg.getParameter<edm::InputTag>("srcRecMuons")),
84  srcRecTracks_(cfg.getParameter<edm::InputTag>("srcRecTracks")),
85  srcCaloTowers_(cfg.getParameter<edm::InputTag>("srcCaloTowers")),
86  srcRecPFCandidates_(cfg.getParameter<edm::InputTag>("srcRecPFCandidates")),
87  srcRecJets_(cfg.getParameter<edm::InputTag>("srcRecJets")),
88  srcTheRecVertex_(cfg.getParameter<edm::InputTag>("srcTheRecVertex")),
89  srcRecVertices_(cfg.getParameter<edm::InputTag>("srcRecVertices")),
90  srcRecVerticesWithBS_(cfg.getParameter<edm::InputTag>("srcRecVerticesWithBS")),
91  srcBeamSpot_(cfg.getParameter<edm::InputTag>("srcBeamSpot")),
92  srcGenDiTaus_(cfg.getParameter<edm::InputTag>("srcGenDiTaus")),
93  dRminSeparation_(cfg.getParameter<double>("dRminSeparation")),
94  srcGenLeg1_(cfg.getParameter<edm::InputTag>("srcGenLeg1")),
95  srcRecLeg1_(cfg.getParameter<edm::InputTag>("srcRecLeg1")),
96  srcGenLeg2_(cfg.getParameter<edm::InputTag>("srcGenLeg2")),
97  srcRecLeg2_(cfg.getParameter<edm::InputTag>("srcRecLeg2")),
98  srcGenParticles_(cfg.getParameter<edm::InputTag>("srcGenParticles")),
99  srcL1ETM_(cfg.getParameter<edm::InputTag>("srcL1ETM")),
100  srcGenCaloMEt_(cfg.getParameter<edm::InputTag>("srcGenCaloMEt")),
101  srcGenPFMEt_(cfg.getParameter<edm::InputTag>("srcGenPFMEt")),
102  srcRecCaloMEt_(cfg.getParameter<edm::InputTag>("srcRecCaloMEt")),
103  srcRecPFMEt_(cfg.getParameter<edm::InputTag>("srcRecPFMEt")),
104  srcMuonsBeforeRad_(cfg.getParameter<edm::InputTag>("srcMuonsBeforeRad")),
105  srcMuonsAfterRad_(cfg.getParameter<edm::InputTag>("srcMuonsAfterRad")),
106  srcMuonRadCorrWeight_(cfg.getParameter<edm::InputTag>("srcMuonRadCorrWeight")),
107  srcMuonRadCorrWeightUp_(cfg.getParameter<edm::InputTag>("srcMuonRadCorrWeightUp")),
108  srcMuonRadCorrWeightDown_(cfg.getParameter<edm::InputTag>("srcMuonRadCorrWeightDown")),
109  srcOtherWeights_(cfg.getParameter<vInputTag>("srcOtherWeights")),
110  srcGenFilterInfo_(cfg.getParameter<edm::InputTag>("srcGenFilterInfo")),
111  dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
112  replacedMuonPtThresholdHigh_(cfg.getParameter<double>("replacedMuonPtThresholdHigh")),
113  replacedMuonPtThresholdLow_(cfg.getParameter<double>("replacedMuonPtThresholdLow"))
114 {
115  typedef std::pair<int, int> pint;
116  std::vector<pint> jetBins;
117  jetBins.push_back(pint(-1, -1));
118  jetBins.push_back(pint(0, 0));
119  jetBins.push_back(pint(1, 1));
120  jetBins.push_back(pint(2, 2));
121  jetBins.push_back(pint(3, 1000));
122  for ( std::vector<pint>::const_iterator jetBin = jetBins.begin();
123  jetBin != jetBins.end(); ++jetBin ) {
124 //--- setup electron Pt, eta and phi distributions;
125 // electron id & isolation and trigger efficiencies
126  setupLeptonDistribution(jetBin->first, jetBin->second, cfg, "electronDistributions", electronDistributions_);
127  setupElectronDistributionExtra(jetBin->first, jetBin->second, cfg, "electronDistributions", electronDistributionsExtra_);
128  setupLeptonEfficiency(jetBin->first, jetBin->second, cfg, "electronEfficiencies", electronEfficiencies_);
129  setupLeptonEfficiency(jetBin->first, jetBin->second, cfg, "gsfElectronEfficiencies", gsfElectronEfficiencies_);
130  setupLeptonL1TriggerEfficiency(jetBin->first, jetBin->second, cfg, "electronL1TriggerEfficiencies", electronL1TriggerEfficiencies_);
131 
132 //--- setup muon Pt, eta and phi distributions;
133 // muon id & isolation and trigger efficiencies
134  setupLeptonDistribution(jetBin->first, jetBin->second, cfg, "muonDistributions", muonDistributions_);
135  setupLeptonEfficiency(jetBin->first, jetBin->second, cfg, "muonEfficiencies", muonEfficiencies_);
136  setupLeptonL1TriggerEfficiency(jetBin->first, jetBin->second, cfg, "muonL1TriggerEfficiencies", muonL1TriggerEfficiencies_);
137 
138 //--- setup tau Pt, eta and phi distributions;
139 // tau id efficiency
140  setupLeptonDistribution(jetBin->first, jetBin->second, cfg, "tauDistributions", tauDistributions_);
141  setupTauDistributionExtra(jetBin->first, jetBin->second, cfg, "tauDistributions", tauDistributionsExtra_);
142  setupLeptonEfficiency(jetBin->first, jetBin->second, cfg, "tauEfficiencies", tauEfficiencies_);
143 
144 //--- setup Pt, eta and phi distributions of L1Extra objects
145 // (electrons, muons, tau-jet, central and forward jets)
146  setupL1ExtraObjectDistribution(jetBin->first, jetBin->second, cfg, "l1ElectronDistributions", l1ElectronDistributions_);
147  setupL1ExtraObjectDistribution(jetBin->first, jetBin->second, cfg, "l1MuonDistributions", l1MuonDistributions_);
148  setupL1ExtraObjectDistribution(jetBin->first, jetBin->second, cfg, "l1TauDistributions", l1TauDistributions_);
149  setupL1ExtraObjectDistribution(jetBin->first, jetBin->second, cfg, "l1CentralJetDistributions", l1CentralJetDistributions_);
150  setupL1ExtraObjectDistribution(jetBin->first, jetBin->second, cfg, "l1ForwardJetDistributions", l1ForwardJetDistributions_);
151 
152 //--- setup MET Pt and phi distributions;
153 // efficiency of L1 (Calo)MET trigger requirement
154  setupMEtDistribution(jetBin->first, jetBin->second, cfg, "metDistributions", metDistributions_);
155  setupMEtL1TriggerEfficiency(jetBin->first, jetBin->second, cfg, "metL1TriggerEfficiencies", metL1TriggerEfficiencies_);
156  }
157 
158  verbosity_ = ( cfg.exists("verbosity") ) ?
159  cfg.getParameter<int>("verbosity") : 0;
160 }
161 
163 {
164  for ( std::vector<plotEntryTypeEvtWeight*>::iterator it = evtWeightPlotEntries_.begin();
165  it != evtWeightPlotEntries_.end(); ++it ) {
166  delete (*it);
167  }
168 
169  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator it = muonRadCorrUncertaintyPlotEntries_beforeRad_.begin();
171  delete (*it);
172  }
173  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator it = muonRadCorrUncertaintyPlotEntries_afterRad_.begin();
174  it != muonRadCorrUncertaintyPlotEntries_afterRad_.end(); ++it ) {
175  delete (*it);
176  }
177  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator it = muonRadCorrUncertaintyPlotEntries_afterRadAndCorr_.begin();
179  delete (*it);
180  }
181 
182  for ( std::vector<plotEntryTypeL1ETM*>::iterator it = l1ETMplotEntries_.begin();
183  it != l1ETMplotEntries_.end(); ++it ) {
184  delete (*it);
185  }
186 
205 }
206 
207 template <typename T>
209  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<leptonDistributionT<T>*>& leptonDistributions)
210 {
211  if ( cfg.exists(keyword) ) {
212  edm::VParameterSet cfgLeptonDistributions = cfg.getParameter<edm::VParameterSet>(keyword);
213  for ( edm::VParameterSet::const_iterator cfgLeptonDistribution = cfgLeptonDistributions.begin();
214  cfgLeptonDistribution != cfgLeptonDistributions.end(); ++cfgLeptonDistribution ) {
215  edm::InputTag srcGen = cfgLeptonDistribution->getParameter<edm::InputTag>("srcGen");
216  std::string cutGen = cfgLeptonDistribution->exists("cutGen") ?
217  cfgLeptonDistribution->getParameter<std::string>("cutGen") : "";
218  edm::InputTag srcRec = cfgLeptonDistribution->getParameter<edm::InputTag>("srcRec");
219  std::string cutRec = cfgLeptonDistribution->exists("cutRec") ?
220  cfgLeptonDistribution->getParameter<std::string>("cutRec") : "";
221  double dRmatch = cfgLeptonDistribution->exists("dRmatch") ?
222  cfgLeptonDistribution->getParameter<double>("dRmatch") : 0.3;
223  std::string dqmDirectory = dqmDirectory_full(cfgLeptonDistribution->getParameter<std::string>("dqmDirectory"));
224  leptonDistributionT<T>* leptonDistribution = new leptonDistributionT<T>(minJets, maxJets, srcGen, cutGen, srcRec, cutRec, dRmatch, dqmDirectory);
225  leptonDistributions.push_back(leptonDistribution);
226  }
227  }
228 }
229 
231  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<electronDistributionExtra*>& electronDistributionsExtra)
232 {
233  if ( cfg.exists(keyword) ) {
234  edm::VParameterSet cfgLeptonDistributions = cfg.getParameter<edm::VParameterSet>(keyword);
235  for ( edm::VParameterSet::const_iterator cfgLeptonDistribution = cfgLeptonDistributions.begin();
236  cfgLeptonDistribution != cfgLeptonDistributions.end(); ++cfgLeptonDistribution ) {
237  edm::InputTag srcGen = cfgLeptonDistribution->getParameter<edm::InputTag>("srcGen");
238  std::string cutGen = cfgLeptonDistribution->exists("cutGen") ?
239  cfgLeptonDistribution->getParameter<std::string>("cutGen") : "";
240  edm::InputTag srcRec = cfgLeptonDistribution->getParameter<edm::InputTag>("srcRec");
241  std::string cutRec = cfgLeptonDistribution->exists("cutRec") ?
242  cfgLeptonDistribution->getParameter<std::string>("cutRec") : "";
243  double dRmatch = cfgLeptonDistribution->exists("dRmatch") ?
244  cfgLeptonDistribution->getParameter<double>("dRmatch") : 0.3;
245  std::string dqmDirectory = dqmDirectory_full(cfgLeptonDistribution->getParameter<std::string>("dqmDirectory"));
246  electronDistributionExtra* electronDistribution = new electronDistributionExtra(minJets, maxJets, srcGen, cutGen, srcRec, cutRec, dRmatch, dqmDirectory, srcTheRecVertex_);
247  electronDistributionsExtra.push_back(electronDistribution);
248  }
249  }
250 }
251 
253  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<tauDistributionExtra*>& tauDistributionsExtra)
254 {
255  if ( cfg.exists(keyword) ) {
256  edm::VParameterSet cfgLeptonDistributions = cfg.getParameter<edm::VParameterSet>(keyword);
257  for ( edm::VParameterSet::const_iterator cfgLeptonDistribution = cfgLeptonDistributions.begin();
258  cfgLeptonDistribution != cfgLeptonDistributions.end(); ++cfgLeptonDistribution ) {
259  edm::InputTag srcGen = cfgLeptonDistribution->getParameter<edm::InputTag>("srcGen");
260  std::string cutGen = cfgLeptonDistribution->exists("cutGen") ?
261  cfgLeptonDistribution->getParameter<std::string>("cutGen") : "";
262  edm::InputTag srcRec = cfgLeptonDistribution->getParameter<edm::InputTag>("srcRec");
263  std::string cutRec = cfgLeptonDistribution->exists("cutRec") ?
264  cfgLeptonDistribution->getParameter<std::string>("cutRec") : "";
265  double dRmatch = cfgLeptonDistribution->exists("dRmatch") ?
266  cfgLeptonDistribution->getParameter<double>("dRmatch") : 0.3;
267  std::string dqmDirectory = dqmDirectory_full(cfgLeptonDistribution->getParameter<std::string>("dqmDirectory"));
268  tauDistributionExtra* tauDistribution = new tauDistributionExtra(minJets, maxJets, srcGen, cutGen, srcRec, cutRec, dRmatch, dqmDirectory);
269  tauDistributionsExtra.push_back(tauDistribution);
270  }
271  }
272 }
273 
274 template <typename T>
276  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<leptonEfficiencyT<T>*>& leptonEfficiencies)
277 {
278  if ( cfg.exists(keyword) ) {
279  edm::VParameterSet cfgLeptonEfficiencies = cfg.getParameter<edm::VParameterSet>(keyword);
280  for ( edm::VParameterSet::const_iterator cfgLeptonEfficiency = cfgLeptonEfficiencies.begin();
281  cfgLeptonEfficiency != cfgLeptonEfficiencies.end(); ++cfgLeptonEfficiency ) {
282  edm::InputTag srcGen = cfgLeptonEfficiency->getParameter<edm::InputTag>("srcGen");
283  std::string cutGen = cfgLeptonEfficiency->exists("cutGen") ?
284  cfgLeptonEfficiency->getParameter<std::string>("cutGen") : "";
285  edm::InputTag srcRec = cfgLeptonEfficiency->getParameter<edm::InputTag>("srcRec");
286  std::string cutRec = cfgLeptonEfficiency->exists("cutRec") ?
287  cfgLeptonEfficiency->getParameter<std::string>("cutRec") : "";
288  double dRmatch = cfgLeptonEfficiency->exists("dRmatch") ?
289  cfgLeptonEfficiency->getParameter<double>("dRmatch") : 0.3;
290  std::string dqmDirectory = dqmDirectory_full(cfgLeptonEfficiency->getParameter<std::string>("dqmDirectory"));
291  leptonEfficiencyT<T>* leptonEfficiency = new leptonEfficiencyT<T>(minJets, maxJets, srcGen, cutGen, srcRec, cutRec, dRmatch, dqmDirectory);
292  leptonEfficiencies.push_back(leptonEfficiency);
293  }
294  }
295 }
296 
297 template <typename T1, typename T2>
299  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<leptonL1TriggerEfficiencyT1T2<T1,T2>*>& leptonL1TriggerEfficiencies)
300 {
301  if ( cfg.exists(keyword) ) {
302  edm::VParameterSet cfgLeptonL1TriggerEfficiencies = cfg.getParameter<edm::VParameterSet>(keyword);
303  for ( edm::VParameterSet::const_iterator cfgLeptonL1TriggerEfficiency = cfgLeptonL1TriggerEfficiencies.begin();
304  cfgLeptonL1TriggerEfficiency != cfgLeptonL1TriggerEfficiencies.end(); ++cfgLeptonL1TriggerEfficiency ) {
305  edm::InputTag srcRef = cfgLeptonL1TriggerEfficiency->getParameter<edm::InputTag>("srcRef");
306  std::string cutRef = cfgLeptonL1TriggerEfficiency->exists("cutRef") ?
307  cfgLeptonL1TriggerEfficiency->getParameter<std::string>("cutRef") : "";
308  edm::InputTag srcL1 = cfgLeptonL1TriggerEfficiency->getParameter<edm::InputTag>("srcL1");
309  std::string cutL1 = cfgLeptonL1TriggerEfficiency->getParameter<std::string>("cutL1");
310  double dRmatch = cfgLeptonL1TriggerEfficiency->exists("dRmatch") ?
311  cfgLeptonL1TriggerEfficiency->getParameter<double>("dRmatch") : 0.3;
312  std::string dqmDirectory = dqmDirectory_full(cfgLeptonL1TriggerEfficiency->getParameter<std::string>("dqmDirectory"));
313  leptonL1TriggerEfficiencyT1T2<T1,T2>* leptonL1TriggerEfficiency = new leptonL1TriggerEfficiencyT1T2<T1,T2>(minJets, maxJets, srcRef, cutRef, srcL1, cutL1, dRmatch, dqmDirectory);
314  leptonL1TriggerEfficiencies.push_back(leptonL1TriggerEfficiency);
315  }
316  }
317 }
318 
319 template <typename T>
321  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<l1ExtraObjectDistributionT<T>*>& l1ExtraObjectDistributions)
322 {
323  if ( cfg.exists(keyword) ) {
324  edm::VParameterSet cfgL1ExtraObjectDistributions = cfg.getParameter<edm::VParameterSet>(keyword);
325  for ( edm::VParameterSet::const_iterator cfgL1ExtraObjectDistribution = cfgL1ExtraObjectDistributions.begin();
326  cfgL1ExtraObjectDistribution != cfgL1ExtraObjectDistributions.end(); ++cfgL1ExtraObjectDistribution ) {
327  edm::InputTag src = cfgL1ExtraObjectDistribution->getParameter<edm::InputTag>("src");
328  std::string cut = cfgL1ExtraObjectDistribution->exists("cut") ?
329  cfgL1ExtraObjectDistribution->getParameter<std::string>("cut") : "";
330  std::string dqmDirectory = dqmDirectory_full(cfgL1ExtraObjectDistribution->getParameter<std::string>("dqmDirectory"));
331  l1ExtraObjectDistributionT<T>* l1ExtraObjectDistribution = new l1ExtraObjectDistributionT<T>(minJets, maxJets, src, cut, dqmDirectory);
332  l1ExtraObjectDistributions.push_back(l1ExtraObjectDistribution);
333  }
334  }
335 }
336 
338  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<metDistributionType*>& metDistributions)
339 {
340  if ( cfg.exists(keyword) ) {
341  edm::VParameterSet cfgMEtDistributions = cfg.getParameter<edm::VParameterSet>(keyword);
342  for ( edm::VParameterSet::const_iterator cfgMEtDistribution = cfgMEtDistributions.begin();
343  cfgMEtDistribution != cfgMEtDistributions.end(); ++cfgMEtDistribution ) {
344  edm::InputTag srcGen = cfgMEtDistribution->getParameter<edm::InputTag>("srcGen");
345  edm::InputTag srcRec = cfgMEtDistribution->getParameter<edm::InputTag>("srcRec");
346  edm::InputTag srcGenZs = cfgMEtDistribution->getParameter<edm::InputTag>("srcGenZs");
347  std::string dqmDirectory = dqmDirectory_full(cfgMEtDistribution->getParameter<std::string>("dqmDirectory"));
348  metDistributionType* metDistribution = new metDistributionType(minJets, maxJets, srcGen, srcRec, srcGenZs, dqmDirectory);
349  metDistributions.push_back(metDistribution);
350  }
351  }
352 }
353 
355  const edm::ParameterSet& cfg, const std::string& keyword, std::vector<metL1TriggerEfficiencyType*>& metL1TriggerEfficiencies)
356 {
357  if ( cfg.exists(keyword) ) {
358  edm::VParameterSet cfgMEtL1TriggerEfficiencies = cfg.getParameter<edm::VParameterSet>(keyword);
359  for ( edm::VParameterSet::const_iterator cfgMEtL1TriggerEfficiency = cfgMEtL1TriggerEfficiencies.begin();
360  cfgMEtL1TriggerEfficiency != cfgMEtL1TriggerEfficiencies.end(); ++cfgMEtL1TriggerEfficiency ) {
361  edm::InputTag srcRef = cfgMEtL1TriggerEfficiency->getParameter<edm::InputTag>("srcRef");
362  edm::InputTag srcL1 = cfgMEtL1TriggerEfficiency->getParameter<edm::InputTag>("srcL1");
363  double cutL1Et = cfgMEtL1TriggerEfficiency->getParameter<double>("cutL1Et");
364  double cutL1Pt = cfgMEtL1TriggerEfficiency->getParameter<double>("cutL1Pt");
365  std::string dqmDirectory = dqmDirectory_full(cfgMEtL1TriggerEfficiency->getParameter<std::string>("dqmDirectory"));
366  metL1TriggerEfficiencyType* metL1TriggerEfficiency = new metL1TriggerEfficiencyType(minJets, maxJets, srcRef, srcL1, cutL1Et, cutL1Pt, dqmDirectory);
367  metL1TriggerEfficiencies.push_back(metL1TriggerEfficiency);
368  }
369  }
370 }
371 
373 {
374  if ( !edm::Service<DQMStore>().isAvailable() )
375  throw cms::Exception("MuonRadiationAnalyzer::beginJob")
376  << "Failed to access dqmStore !!\n";
377 
378  DQMStore& dqmStore = (*edm::Service<DQMStore>());
379  dqmStore.setCurrentFolder(dqmDirectory_.data());
380 
381 //--- book all histograms
382  histogramEventCounter_ = dqmStore.book1D("EventCounter", "EventCounter", 1, -0.5, 1.5);
383 
384  histogramGenFilterEfficiency_ = dqmStore.book1D("genFilterEfficiency", "genFilterEfficiency", 102, -0.01, 1.01);
385 
386  histogramRotationAngleMatrix_ = dqmStore.book2D("rfRotationAngleMatrix", "rfRotationAngleMatrix", 2, -0.5, 1.5, 2, -0.5, 1.5);
387  histogramRotationLegPlusDeltaR_ = dqmStore.book1D("rfRotationLegPlusDeltaR", "rfRotationLegPlusDeltaR", 101, -0.05, 10.05);
388  histogramRotationLegMinusDeltaR_ = dqmStore.book1D("rfRotationLegMinusDeltaR", "rfRotationLegMinusDeltaR", 101, -0.05, 10.05);
389  histogramPhiRotLegPlus_ = dqmStore.book1D("rfPhiRotLegPlus", "rfPhiRotLegPlus", 72, -TMath::Pi(), +TMath::Pi());
390  histogramPhiRotLegMinus_ = dqmStore.book1D("rfPhiRotLegMinus", "rfPhiRotLegMinus", 72, -TMath::Pi(), +TMath::Pi());
391 
392  histogramNumTracksPtGt5_ = dqmStore.book1D("numTracksPtGt5", "numTracksPtGt5", 50, -0.5, 49.5);
393  histogramNumTracksPtGt10_ = dqmStore.book1D("numTracksPtGt10", "numTracksPtGt10", 50, -0.5, 49.5);
394  histogramNumTracksPtGt20_ = dqmStore.book1D("numTracksPtGt20", "numTracksPtGt20", 50, -0.5, 49.5);
395  histogramNumTracksPtGt30_ = dqmStore.book1D("numTracksPtGt30", "numTracksPtGt30", 50, -0.5, 49.5);
396  histogramNumTracksPtGt40_ = dqmStore.book1D("numTracksPtGt40", "numTracksPtGt40", 50, -0.5, 49.5);
397 
398  histogramNumGlobalMuons_ = dqmStore.book1D("numGlobalMuons", "numGlobalMuons", 20, -0.5, 19.5);
399  histogramNumStandAloneMuons_ = dqmStore.book1D("numStandAloneMuons", "numStandAloneMuons", 20, -0.5, 19.5);
400  histogramNumPFMuons_ = dqmStore.book1D("numPFMuons", "numPFMuons", 20, -0.5, 19.5);
401 
402  histogramNumChargedPFCandsPtGt5_ = dqmStore.book1D("numChargedPFCandsPtGt5", "numChargedPFCandsPtGt5", 50, -0.5, 49.5);
403  histogramNumChargedPFCandsPtGt10_ = dqmStore.book1D("numChargedPFCandsPtGt10", "numChargedPFCandsPtGt10", 50, -0.5, 49.5);
404  histogramNumChargedPFCandsPtGt20_ = dqmStore.book1D("numChargedPFCandsPtGt20", "numChargedPFCandsPtGt20", 50, -0.5, 49.5);
405  histogramNumChargedPFCandsPtGt30_ = dqmStore.book1D("numChargedPFCandsPtGt30", "numChargedPFCandsPtGt30", 50, -0.5, 49.5);
406  histogramNumChargedPFCandsPtGt40_ = dqmStore.book1D("numChargedPFCandsPtGt40", "numChargedPFCandsPtGt40", 50, -0.5, 49.5);
407 
408  histogramNumNeutralPFCandsPtGt5_ = dqmStore.book1D("numNeutralPFCandsPtGt5", "numNeutralPFCandsPtGt5", 50, -0.5, 49.5);
409  histogramNumNeutralPFCandsPtGt10_ = dqmStore.book1D("numNeutralPFCandsPtGt10", "numNeutralPFCandsPtGt10", 50, -0.5, 49.5);
410  histogramNumNeutralPFCandsPtGt20_ = dqmStore.book1D("numNeutralPFCandsPtGt20", "numNeutralPFCandsPtGt20", 50, -0.5, 49.5);
411  histogramNumNeutralPFCandsPtGt30_ = dqmStore.book1D("numNeutralPFCandsPtGt30", "numNeutralPFCandsPtGt30", 50, -0.5, 49.5);
412  histogramNumNeutralPFCandsPtGt40_ = dqmStore.book1D("numNeutralPFCandsPtGt40", "numNeutralPFCandsPtGt40", 50, -0.5, 49.5);
413 
414  histogramRawJetPt_ = dqmStore.book1D("rawJetPt", "rawJetPt", 250, 0., 250.);
415  histogramRawJetPtAbsEtaLt2_5_ = dqmStore.book1D("rawJetPtAbsEtaLt2_5", "rawJetPtAbsEtaLt2_5", 250, 0., 250.);
416  histogramRawJetPtAbsEta2_5to4_5_ = dqmStore.book1D("rawJetPtAbsEta2_5to4_5", "rawJetPtAbsEta2_5to4_5", 250, 0., 250.);
417  histogramRawJetEtaPtGt20_ = dqmStore.book1D("rawJetEtaPtGt20", "rawJetEtaPtGt20", 198, -9.9, +9.9);
418  histogramRawJetEtaPtGt30_ = dqmStore.book1D("rawJetEtaPtGt30", "rawJetEtaPtGt30", 198, -9.9, +9.9);
419  histogramNumJetsRawPtGt20_ = dqmStore.book1D("numJetsRawPtGt20", "numJetsRawPtGt20", 50, -0.5, 49.5);
420  histogramNumJetsRawPtGt20AbsEtaLt2_5_ = dqmStore.book1D("numJetsRawPtGt20AbsEtaLt2_5", "numJetsRawPtGt20AbsEtaLt2_5", 50, -0.5, 49.5);
421  histogramNumJetsRawPtGt20AbsEta2_5to4_5_ = dqmStore.book1D("numJetsRawPtGt20AbsEta2_5to4_5", "numJetsRawPtGt20AbsEta2_5to4_5", 50, -0.5, 49.5);
422  histogramNumJetsRawPtGt30_ = dqmStore.book1D("numJetsRawPtGt30", "numJetsRawPtGt30", 50, -0.5, 49.5);
423  histogramNumJetsRawPtGt30AbsEtaLt2_5_ = dqmStore.book1D("numJetsRawPtGt30AbsEtaLt2_5", "numJetsRawPtGt30AbsEtaLt2_5", 50, -0.5, 49.5);
424  histogramNumJetsRawPtGt30AbsEta2_5to4_5_ = dqmStore.book1D("numJetsRawPtGt30AbsEta2_5to4_5", "numJetsRawPtGt30AbsEta2_5to4_5", 50, -0.5, 49.5);
425  histogramCorrJetPt_ = dqmStore.book1D("corrJetPt", "corrJetPt", 250, 0., 250.);
426  histogramCorrJetPtAbsEtaLt2_5_ = dqmStore.book1D("corrJetPtAbsEtaLt2_5", "corrJetPtAbsEtaLt2_5", 250, 0., 250.);
427  histogramCorrJetPtAbsEta2_5to4_5_ = dqmStore.book1D("corrJetPtAbsEta2_5to4_5", "corrJetPtAbsEta2_5to4_5", 250, 0., 250.);
428  histogramCorrJetEtaPtGt20_ = dqmStore.book1D("corrJetEtaPtGt20", "corrJetEtaPtGt20", 198, -9.9, +9.9);
429  histogramCorrJetEtaPtGt30_ = dqmStore.book1D("corrJetEtaPtGt30", "corrJetEtaPtGt30", 198, -9.9, +9.9);
430  histogramNumJetsCorrPtGt20_ = dqmStore.book1D("numJetsCorrPtGt20", "numJetsCorrPtGt20", 20, -0.5, 19.5);
431  histogramNumJetsCorrPtGt20AbsEtaLt2_5_ = dqmStore.book1D("numJetsCorrPtGt20AbsEtaLt2_5", "numJetsCorrPtGt20AbsEtaLt2_5", 20, -0.5, 19.5);
432  histogramNumJetsCorrPtGt20AbsEta2_5to4_5_ = dqmStore.book1D("numJetsCorrPtGt20AbsEta2_5to4_5", "numJetsCorrPtGt20AbsEta2_5to4_5", 20, -0.5, 19.5);
433  histogramNumJetsCorrPtGt30_ = dqmStore.book1D("numJetsCorrPtGt30", "numJetsCorrPtGt30", 20, -0.5, 19.5);
434  histogramNumJetsCorrPtGt30AbsEtaLt2_5_ = dqmStore.book1D("numJetsCorrPtGt30AbsEtaLt2_5", "numJetsCorrPtGt30AbsEtaLt2_5", 20, -0.5, 19.5);
435  histogramNumJetsCorrPtGt30AbsEta2_5to4_5_ = dqmStore.book1D("numJetsCorrPtGt30AbsEta2_5to4_5", "numJetsCorrPtGt30AbsEta2_5to4_5", 20, -0.5, 19.5);
436 
437  histogramTheRecVertexX_ = dqmStore.book1D("theRecVertexX", "theRecVertexX", 2000, -1., +1.);
438  histogramTheRecVertexY_ = dqmStore.book1D("theRecVertexY", "theRecVertexY", 2000, -1., +1.);
439  histogramTheRecVertexZ_ = dqmStore.book1D("theRecVertexZ", "theRecVertexZ", 500, -25., +25.);
440  histogramRecVertexX_ = dqmStore.book1D("recVertexX", "recVertexX", 2000, -1., +1.);
441  histogramRecVertexY_ = dqmStore.book1D("recVertexY", "recVertexY", 2000, -1., +1.);
442  histogramRecVertexZ_ = dqmStore.book1D("recVertexZ", "recVertexZ", 500, -25., +25.);
443  histogramNumRecVertices_ = dqmStore.book1D("numRecVertices", "numRecVertices", 50, -0.5, +49.5);
444  histogramRecVertexWithBSx_ = dqmStore.book1D("recVertexWithBSx", "recVertexWithBSx", 2000, -1., +1.);
445  histogramRecVertexWithBSy_ = dqmStore.book1D("recVertexWithBSy", "recVertexWithBSy", 2000, -1., +1.);
446  histogramRecVertexWithBSz_ = dqmStore.book1D("recVertexWithBSz", "recVertexWithBSz", 500, -25., +25.);
447  histogramNumRecVerticesWithBS_ = dqmStore.book1D("numRecVerticesWithBS", "numRecVerticesWithBS", 50, -0.5, +49.5);
448 
449  histogramBeamSpotX_ = dqmStore.book1D("beamSpotX", "beamSpotX", 2000, -1., +1.);
450  histogramBeamSpotY_ = dqmStore.book1D("beamSpotY", "beamSpotY", 2000, -1., +1.);
451 
452  histogramGenDiTauPt_ = dqmStore.book1D("genDiTauPt", "genDiTauPt", 250, 0., 250.);
453  histogramGenDiTauEta_ = dqmStore.book1D("genDiTauEta", "genDiTauEta", 198, -9.9, +9.9);
454  histogramGenDiTauPhi_ = dqmStore.book1D("genDiTauPhi", "genDiTauPhi", 72, -TMath::Pi(), +TMath::Pi());
455  histogramGenDiTauMass_ = dqmStore.book1D("genDiTauMass", "genDiTauMass", 500, 0., 500.);
456  histogramGenDeltaPhiLeg1Leg2_ = dqmStore.book1D("genDeltaPhiLeg1Leg2", "genDeltaPhiLeg1Leg2", 180, 0., +TMath::Pi());
457  histogramGenDiTauDecayAngle_ = dqmStore.book1D("genDiTauDecayAngle", "genDiTauDecayAngle", 180, 0., +TMath::Pi());
458 
459  histogramGenVisDiTauPt_ = dqmStore.book1D("genVisDiTauPt", "genVisDiTauPt", 250, 0., 250.);
460  histogramGenVisDiTauEta_ = dqmStore.book1D("genVisDiTauEta", "genVisDiTauEta", 198, -9.9, +9.9);
461  histogramGenVisDiTauPhi_ = dqmStore.book1D("genVisDiTauPhi", "genVisDiTauPhi", 72, -TMath::Pi(), +TMath::Pi());
462  histogramGenVisDiTauMass_ = dqmStore.book1D("genVisDiTauMass", "genVisDiTauMass", 500, 0., 500.);
463  histogramGenVisDeltaPhiLeg1Leg2_ = dqmStore.book1D("genVisDeltaPhiLeg1Leg2", "genVisDeltaPhiLeg1Leg2", 180, 0., +TMath::Pi());
464 
465  histogramRecVisDiTauPt_ = dqmStore.book1D("recVisDiTauPt", "recVisDiTauPt", 250, 0., 250.);
466  histogramRecVisDiTauEta_ = dqmStore.book1D("recVisDiTauEta", "recVisDiTauEta", 198, -9.9, +9.9);
467  histogramRecVisDiTauPhi_ = dqmStore.book1D("recVisDiTauPhi", "recVisDiTauPhi", 72, -TMath::Pi(), +TMath::Pi());
468  histogramRecVisDiTauMass_ = dqmStore.book1D("recVisDiTauMass", "recVisDiTauMass", 500, 0., 500.);
469  histogramRecVisDeltaPhiLeg1Leg2_ = dqmStore.book1D("recVisDeltaPhiLeg1Leg2", "recVisDeltaPhiLeg1Leg2", 180, 0., +TMath::Pi());
470 
471  histogramGenTau1Pt_ = dqmStore.book1D("genTau1Pt", "genTau1Pt", 250, 0., 250.);
472  histogramGenTau1Eta_ = dqmStore.book1D("genTau1Eta", "genTau1Eta", 198, -9.9, +9.9);
473  histogramGenTau1Phi_ = dqmStore.book1D("genTau1Phi", "genTau1Phi", 72, -TMath::Pi(), +TMath::Pi());
474  histogramGenLeg1Pt_ = dqmStore.book1D("genLeg1Pt", "genLeg1Pt", 250, 0., 250.);
475  histogramGenLeg1Eta_ = dqmStore.book1D("genLeg1Eta", "genLeg1Eta", 198, -9.9, +9.9);
476  histogramGenLeg1Phi_ = dqmStore.book1D("genLeg1Phi", "genLeg1Phi", 72, -TMath::Pi(), +TMath::Pi());
477  histogramGenLeg1X_ = dqmStore.book1D("genLeg1X", "X_{1}^{gen}", 102, -0.01, 1.01);
478  histogramGenLeg1XforGenLeg2X0_00to0_25_ = dqmStore.book1D("genLeg1XforGenLeg2X0_00to0_25", "X_{1}^{gen} (0.00 < X_{2}^{gen} < 0.25)", 102, -0.01, 1.01);
479  histogramGenLeg1XforGenLeg2X0_25to0_50_ = dqmStore.book1D("genLeg1XforGenLeg2X0_25to0_50", "X_{1}^{gen} (0.25 < X_{2}^{gen} < 0.50)", 102, -0.01, 1.01);
480  histogramGenLeg1XforGenLeg2X0_50to0_75_ = dqmStore.book1D("genLeg1XforGenLeg2X0_50to0_75", "X_{1}^{gen} (0.50 < X_{2}^{gen} < 0.75)", 102, -0.01, 1.01);
481  histogramGenLeg1XforGenLeg2X0_75to1_00_ = dqmStore.book1D("genLeg1XforGenLeg2X0_75to1_00", "X_{1}^{gen} (0.75 < X_{2}^{gen} < 1.00)", 102, -0.01, 1.01);
482  histogramGenLeg1Mt_ = dqmStore.book1D("genLeg1Mt", "genLeg1Mt", 250, 0., 250.);
483  histogramRecLeg1X_ = dqmStore.book1D("recLeg1X", "recLeg1X", 102, -0.01, 1.01);
484  histogramRecLeg1PFMt_ = dqmStore.book1D("recLeg1PFMt", "recLeg1PFMt", 250, 0., 250.);
485  histogramGenTau2Pt_ = dqmStore.book1D("genTau2Pt", "genTau2Pt", 250, 0., 250.);
486  histogramGenTau2Eta_ = dqmStore.book1D("genTau2Eta", "genTau2Eta", 198, -9.9, +9.9);
487  histogramGenTau2Phi_ = dqmStore.book1D("genTau2Phi", "genTau2Phi", 72, -TMath::Pi(), +TMath::Pi());
488  histogramGenLeg2Pt_ = dqmStore.book1D("genLeg2Pt", "genLeg2Pt", 250, 0., 250.);
489  histogramGenLeg2Eta_ = dqmStore.book1D("genLeg2Eta", "genLeg2Eta", 198, -9.9, +9.9);
490  histogramGenLeg2Phi_ = dqmStore.book1D("genLeg2Phi", "genLeg2Phi", 72, -TMath::Pi(), +TMath::Pi());
491  histogramGenLeg2X_ = dqmStore.book1D("genLeg2X", "X_{2}^{gen}", 102, -0.01, 1.01);
492  histogramGenLeg2XforGenLeg1X0_00to0_25_ = dqmStore.book1D("genLeg2XforGenLeg1X0_00to0_25", "X_{2}^{gen} (0.00 < X_{1}^{gen} < 0.25)", 102, -0.01, 1.01);
493  histogramGenLeg2XforGenLeg1X0_25to0_50_ = dqmStore.book1D("genLeg2XforGenLeg1X0_25to0_50", "X_{2}^{gen} (0.25 < X_{1}^{gen} < 0.50)", 102, -0.01, 1.01);
494  histogramGenLeg2XforGenLeg1X0_50to0_75_ = dqmStore.book1D("genLeg2XforGenLeg1X0_50to0_75", "X_{2}^{gen} (0.50 < X_{1}^{gen} < 0.75)", 102, -0.01, 1.01);
495  histogramGenLeg2XforGenLeg1X0_75to1_00_ = dqmStore.book1D("genLeg2XforGenLeg1X0_75to1_00", "X_{2}^{gen} (0.75 < X_{1}^{gen} < 1.00)", 102, -0.01, 1.01);
496  histogramGenLeg2Mt_ = dqmStore.book1D("genLeg2Mt", "genLeg2Mt", 250, 0., 250.);
497  histogramRecLeg2X_ = dqmStore.book1D("recLeg2X", "recLeg2X", 102, -0.01, 1.01);
498  histogramRecLeg2PFMt_ = dqmStore.book1D("recLeg2PFMt", "recLeg2PFMt", 250, 0., 250.);
499 
500  histogramSumGenParticlePt_ = dqmStore.book1D("sumGenParticlePt", "sumGenParticlePt", 250, 0., 250.);
501  histogramSumGenParticlePt_charged_ = dqmStore.book1D("sumGenParticlePt_charged", "sumGenParticlePt_charged", 250, 0., 250.);
502  histogramGenCaloMEt_ = dqmStore.book1D("genCaloMEt", "genCaloMEt", 250, 0., 250.);
503  histogramGenPFMEt_ = dqmStore.book1D("genPFMEt", "genPFMEt", 250, 0., 250.);
504 
505  histogramRecCaloMEtECAL_ = dqmStore.book1D("recCaloMEtECAL", "recCaloMEtECAL", 250, 0., 250.);
506  histogramRecCaloSumEtECAL_ = dqmStore.book1D("recCaloSumEtECAL", "recCaloSumEtECAL", 2500, 0., 2500.);
507  histogramRecCaloMEtHCAL_ = dqmStore.book1D("recCaloMEtHCAL", "recCaloMEtHCAL", 250, 0., 250.);
508  histogramRecCaloSumEtHCAL_ = dqmStore.book1D("recCaloSumEtHCAL", "recCaloSumEtHCAL", 2500, 0., 2500.);
509  histogramRecCaloMEtHF_ = dqmStore.book1D("recCaloMEtHF", "recCaloMEtHF", 250, 0., 250.);
510  histogramRecCaloSumEtHF_ = dqmStore.book1D("recCaloSumEtHF", "recCaloSumEtHF", 2500, 0., 2500.);
511  histogramRecCaloMEtHO_ = dqmStore.book1D("recCaloMEtHO", "recCaloMEtHO", 250, 0., 250.);
512  histogramRecCaloSumEtHO_ = dqmStore.book1D("recCaloSumEtHO", "recCaloSumEtHO", 2500, 0., 2500.);
513 
514  // CV: Record presence in the embedded event of high Pt tracks, PFCandidates and muons
515  // reconstructed near the direction of the replaced muons
516  // (indicating that maybe not all muon signals were removed in the embedding).
517  // Fill product of reconstructed track/PFCandidate/muon charge * charge of replaced muon into histogram.
518  // In case all matches happen just by chance, expect equal number of entries in positive and negative bins
519  histogramWarning_recTrackNearReplacedMuon_ = dqmStore.book1D("Warning_recTrackNearReplacedMuon", "Warning_recTrackNearReplacedMuon", 3, -1.5, +1.5);
520  histogramWarning_recPFCandNearReplacedMuon_ = dqmStore.book1D("Warning_recPFCandNearReplacedMuon", "Warning_recPFCandNearReplacedMuon", 3, -1.5, +1.5);
521  histogramWarning_recMuonNearReplacedMuon_ = dqmStore.book1D("Warning_recMuonNearReplacedMuon", "Warning_recMuonNearReplacedMuon", 3, -1.5, +1.5);
522 
523  for ( vInputTag::const_iterator srcWeight = srcOtherWeights_.begin();
524  srcWeight != srcOtherWeights_.end(); ++srcWeight ) {
525  plotEntryTypeEvtWeight* evtWeightPlotEntry = new plotEntryTypeEvtWeight(*srcWeight, dqmDirectory_);
526  evtWeightPlotEntry->bookHistograms(dqmStore);
527  evtWeightPlotEntries_.push_back(evtWeightPlotEntry);
528  }
529  if ( srcMuonRadCorrWeight_.label() != "" ) {
531  evtWeightPlotEntry->bookHistograms(dqmStore);
532  evtWeightPlotEntries_.push_back(evtWeightPlotEntry);
533  }
534 
535  typedef std::pair<int, int> pint;
536  std::vector<pint> jetBins;
537  jetBins.push_back(pint(-1, -1));
538  jetBins.push_back(pint(0, 0));
539  jetBins.push_back(pint(1, 1));
540  jetBins.push_back(pint(2, 2));
541  jetBins.push_back(pint(3, 1000));
542  for ( std::vector<pint>::const_iterator jetBin = jetBins.begin();
543  jetBin != jetBins.end(); ++jetBin ) {
544  TString dqmDirectory_beforeRad = dqmDirectory_;
545  if ( !dqmDirectory_beforeRad.EndsWith("/") ) dqmDirectory_beforeRad.Append("/");
546  dqmDirectory_beforeRad.Append("beforeRad");
547  plotEntryTypeMuonRadCorrUncertainty* muonRadCorrUncertaintyPlotEntry_beforeRad = new plotEntryTypeMuonRadCorrUncertainty(jetBin->first, jetBin->second, dqmDirectory_beforeRad.Data());
548  muonRadCorrUncertaintyPlotEntry_beforeRad->bookHistograms(dqmStore);
549  muonRadCorrUncertaintyPlotEntries_beforeRad_.push_back(muonRadCorrUncertaintyPlotEntry_beforeRad);
550  TString dqmDirectory_afterRad = dqmDirectory_;
551  if ( !dqmDirectory_afterRad.EndsWith("/") ) dqmDirectory_afterRad.Append("/");
552  dqmDirectory_afterRad.Append("afterRad");
553  plotEntryTypeMuonRadCorrUncertainty* muonRadCorrUncertaintyPlotEntry_afterRad = new plotEntryTypeMuonRadCorrUncertainty(jetBin->first, jetBin->second, dqmDirectory_afterRad.Data());
554  muonRadCorrUncertaintyPlotEntry_afterRad->bookHistograms(dqmStore);
555  muonRadCorrUncertaintyPlotEntries_afterRad_.push_back(muonRadCorrUncertaintyPlotEntry_afterRad);
556  TString dqmDirectory_afterRadAndCorr = dqmDirectory_;
557  if ( !dqmDirectory_afterRadAndCorr.EndsWith("/") ) dqmDirectory_afterRadAndCorr.Append("/");
558  dqmDirectory_afterRadAndCorr.Append("afterRadAndCorr");
559  plotEntryTypeMuonRadCorrUncertainty* muonRadCorrUncertaintyPlotEntry_afterRadAndCorr = new plotEntryTypeMuonRadCorrUncertainty(jetBin->first, jetBin->second, dqmDirectory_afterRadAndCorr.Data());
560  muonRadCorrUncertaintyPlotEntry_afterRadAndCorr->bookHistograms(dqmStore);
561  muonRadCorrUncertaintyPlotEntries_afterRadAndCorr_.push_back(muonRadCorrUncertaintyPlotEntry_afterRadAndCorr);
562  }
565 
566  std::vector<std::string> genTauDecayModes;
567  genTauDecayModes.push_back(std::string("")); // all tau decay modes
568  genTauDecayModes.push_back(std::string("oneProng0Pi0"));
569  genTauDecayModes.push_back(std::string("oneProng1Pi0"));
570  genTauDecayModes.push_back(std::string("oneProng2Pi0"));
571  genTauDecayModes.push_back(std::string("threeProng0Pi0"));
572  genTauDecayModes.push_back(std::string("threeProng1Pi0"));
573  for ( std::vector<std::string>::const_iterator genTauDecayMode = genTauDecayModes.begin();
574  genTauDecayMode != genTauDecayModes.end(); ++genTauDecayMode ) {
576  l1ETMplotEntry->bookHistograms(dqmStore);
577  l1ETMplotEntries_.push_back(l1ETMplotEntry);
578  }
579 
590  bookHistograms(tauEfficiencies_, dqmStore);
598 }
599 
600 namespace
601 {
602  void fillVisPtEtaPhiMassDistributions(const edm::Event& evt,
603  const edm::InputTag& srcLeg1, const edm::InputTag& srcLeg2,
604  MonitorElement* histogram_visDiTauPt, MonitorElement* histogram_visDiTauEta, MonitorElement* histogram_visDiTauPhi, MonitorElement* histogram_visDiTauMass,
605  double evtWeight)
606  {
607  //std::cout << "<fillVisPtEtaPhiMassDistributions>:" << std::endl;
608  //std::cout << " srcLeg1 = " << srcLeg1 << std::endl;
609  //std::cout << " srcLeg2 = " << srcLeg2 << std::endl;
611  edm::Handle<CandidateView> visDecayProducts1;
612  evt.getByLabel(srcLeg1, visDecayProducts1);
613  edm::Handle<CandidateView> visDecayProducts2;
614  evt.getByLabel(srcLeg2, visDecayProducts2);
615  for ( edm::View<reco::Candidate>::const_iterator visDecayProduct1 = visDecayProducts1->begin();
616  visDecayProduct1 != visDecayProducts1->end(); ++visDecayProduct1 ) {
617  for ( edm::View<reco::Candidate>::const_iterator visDecayProduct2 = visDecayProducts2->begin();
618  visDecayProduct2 != visDecayProducts2->end(); ++visDecayProduct2 ) {
619  double dR = deltaR(visDecayProduct1->p4(), visDecayProduct2->p4());
620  if ( dR > 0.3 ) { // protection in case srcLeg1 and srcLeg2 refer to same collection (e.g. both hadronic tau decays)
621  //std::cout << "leg1: Pt = " << visDecayProduct1->pt() << ", phi = " << visDecayProduct1->phi() << " (Px = " << visDecayProduct1->px() << ", Py = " << visDecayProduct1->py() << ")" << std::endl;
622  //std::cout << "leg2: Pt = " << visDecayProduct2->pt() << ", phi = " << visDecayProduct2->phi() << " (Px = " << visDecayProduct2->px() << ", Py = " << visDecayProduct2->py() << ")" << std::endl;
623  reco::Candidate::LorentzVector visDiTauP4 = visDecayProduct1->p4() + visDecayProduct2->p4();
624  histogram_visDiTauPt->Fill(visDiTauP4.pt(), evtWeight);
625  histogram_visDiTauEta->Fill(visDiTauP4.eta(), evtWeight);
626  histogram_visDiTauPhi->Fill(visDiTauP4.phi(), evtWeight);
627  histogram_visDiTauMass->Fill(visDiTauP4.mass(), evtWeight);
628  }
629  }
630  }
631  }
632 
633  double compMt(const reco::Candidate::LorentzVector& visP4, const reco::Candidate::LorentzVector& metP4)
634  {
635  double sumPx = visP4.px() + metP4.px();
636  double sumPy = visP4.py() + metP4.py();
637  double sumEt = TMath::Max(visP4.Et(), visP4.pt()) + TMath::Sqrt(metP4.pt());
638  double mt2 = sumEt*sumEt - (sumPx*sumPx + sumPy*sumPy);
639  if ( mt2 < 0 ) mt2 = 0.;
640  return TMath::Sqrt(mt2);
641  }
642 
643  void fillX1andX2Distributions(const edm::Event& evt,
644  const edm::InputTag& srcGenDiTau, const edm::InputTag& srcLeg1, const edm::InputTag& srcLeg2, const reco::Candidate::LorentzVector& metP4,
645  MonitorElement* histogram_tau1Pt, MonitorElement* histogram_tau1Eta, MonitorElement* histogram_tau1Phi,
646  MonitorElement* histogram_leg1Pt, MonitorElement* histogram_leg1Eta, MonitorElement* histogram_leg1Phi,
647  MonitorElement* histogram_leg1X,
648  MonitorElement* histogram_leg1XforLeg2X0_00to0_25,
649  MonitorElement* histogram_leg1XforLeg2X0_25to0_50,
650  MonitorElement* histogram_leg1XforLeg2X0_50to0_75,
651  MonitorElement* histogram_leg1XforLeg2X0_75to1_00,
652  MonitorElement* histogram_leg1Mt,
653  MonitorElement* histogram_tau2Pt, MonitorElement* histogram_tau2Eta, MonitorElement* histogram_tau2Phi,
654  MonitorElement* histogram_leg2Pt, MonitorElement* histogram_leg2Eta, MonitorElement* histogram_leg2Phi,
655  MonitorElement* histogram_leg2X,
656  MonitorElement* histogram_leg2XforLeg1X0_00to0_25,
657  MonitorElement* histogram_leg2XforLeg1X0_25to0_50,
658  MonitorElement* histogram_leg2XforLeg1X0_50to0_75,
659  MonitorElement* histogram_leg2XforLeg1X0_75to1_00,
660  MonitorElement* histogram_leg2Mt,
661  MonitorElement* histogram_VisDeltaPhiLeg1Leg2,
662  double evtWeight)
663  {
664  //std::cout << "<fillX1andX2Distributions>:" << std::endl;
665  //std::cout << " srcLeg1 = " << srcLeg1.label() << std::endl;
666  //std::cout << " srcLeg2 = " << srcLeg2.label() << std::endl;
668  edm::Handle<CandidateView> genDiTaus;
669  evt.getByLabel(srcGenDiTau, genDiTaus);
670  edm::Handle<CandidateView> visDecayProducts1;
671  evt.getByLabel(srcLeg1, visDecayProducts1);
672  //std::cout << "#visDecayProducts1 = " << visDecayProducts1->size() << std::endl;
673  edm::Handle<CandidateView> visDecayProducts2;
674  evt.getByLabel(srcLeg2, visDecayProducts2);
675  //std::cout << "#visDecayProducts2 = " << visDecayProducts2->size() << std::endl;
676  for ( CandidateView::const_iterator genDiTau = genDiTaus->begin();
677  genDiTau != genDiTaus->end(); ++genDiTau ) {
678  const reco::CompositeCandidate* genDiTau_composite = dynamic_cast<const reco::CompositeCandidate*>(&(*genDiTau));
679  if ( !(genDiTau_composite && genDiTau_composite->numberOfDaughters() == 2) ) continue;
680  const reco::Candidate* genLeg1 = genDiTau_composite->daughter(0);
681  const reco::Candidate* genLeg2 = genDiTau_composite->daughter(1);
682  //std::cout << "genLeg1: Pt = " << genLeg1->pt() << ", eta = " << genLeg1->eta() << ", phi = " << genLeg1->phi() << std::endl;
683  //std::cout << "genLeg2: Pt = " << genLeg2->pt() << ", eta = " << genLeg2->eta() << ", phi = " << genLeg2->phi() << std::endl;
684  //std::cout << "genLeg1+2: mass = " << (genLeg1->p4() + genLeg2->p4()).mass() << std::endl;
685  if ( !(genLeg1 && genLeg1->energy() > 0. &&
686  genLeg2 && genLeg2->energy() > 0.) ) continue;
687  for ( edm::View<reco::Candidate>::const_iterator visDecayProduct1 = visDecayProducts1->begin();
688  visDecayProduct1 != visDecayProducts1->end(); ++visDecayProduct1 ) {
689  double dR1matchedTo1 = deltaR(visDecayProduct1->p4(), genLeg1->p4());
690  double dR1matchedTo2 = deltaR(visDecayProduct1->p4(), genLeg2->p4());
691  for ( edm::View<reco::Candidate>::const_iterator visDecayProduct2 = visDecayProducts2->begin();
692  visDecayProduct2 != visDecayProducts2->end(); ++visDecayProduct2 ) {
693  double dR2matchedTo1 = deltaR(visDecayProduct2->p4(), genLeg1->p4());
694  double dR2matchedTo2 = deltaR(visDecayProduct2->p4(), genLeg2->p4());
696  double X1 = 0.;
698  double X2 = 0.;
699  bool matched = false;
700  if ( dR1matchedTo1 < 0.3 && dR2matchedTo1 > 0.5 &&
701  dR2matchedTo2 < 0.3 && dR1matchedTo2 > 0.5 ) {
702  leg1P4 = genLeg1->p4();
703  X1 = visDecayProduct1->energy()/genLeg1->energy();
704  leg2P4 = genLeg2->p4();
705  X2 = visDecayProduct2->energy()/genLeg2->energy();
706  matched = true;
707  } else if ( dR1matchedTo2 < 0.3 && dR2matchedTo2 > 0.5 &&
708  dR2matchedTo1 < 0.3 && dR1matchedTo1 > 0.5 ) {
709  leg1P4 = genLeg2->p4();
710  X1 = visDecayProduct1->energy()/genLeg2->energy();
711  leg2P4 = genLeg1->p4();
712  X2 = visDecayProduct2->energy()/genLeg1->energy();
713  matched = true;
714  }
715  if ( matched ) {
716  //std::cout << "X1 = " << X1 << ", X2 = " << X2 << std::endl;
717  if ( histogram_tau1Pt && histogram_tau1Eta && histogram_tau1Phi ) {
718  histogram_tau1Pt->Fill(genLeg1->pt(), evtWeight);
719  histogram_tau1Eta->Fill(genLeg1->eta(), evtWeight);
720  histogram_tau1Phi->Fill(genLeg1->phi(), evtWeight);
721  }
722  if ( histogram_leg1Pt && histogram_leg1Eta && histogram_leg1Phi ) {
723  histogram_leg1Pt->Fill(visDecayProduct1->pt(), evtWeight);
724  histogram_leg1Eta->Fill(visDecayProduct1->eta(), evtWeight);
725  histogram_leg1Phi->Fill(visDecayProduct1->phi(), evtWeight);
726  }
727  histogram_leg1X->Fill(X1, evtWeight);
728  if ( histogram_leg1XforLeg2X0_00to0_25 && histogram_leg1XforLeg2X0_25to0_50 && histogram_leg1XforLeg2X0_50to0_75 && histogram_leg1XforLeg2X0_75to1_00 ) {
729  if ( X2 < 0.25 ) histogram_leg1XforLeg2X0_00to0_25->Fill(X1, evtWeight);
730  else if ( X2 < 0.50 ) histogram_leg1XforLeg2X0_25to0_50->Fill(X1, evtWeight);
731  else if ( X2 < 0.75 ) histogram_leg1XforLeg2X0_50to0_75->Fill(X1, evtWeight);
732  else histogram_leg1XforLeg2X0_75to1_00->Fill(X1, evtWeight);
733  }
734  if ( histogram_leg1Mt ) histogram_leg1Mt->Fill(compMt(visDecayProduct1->p4(), metP4), evtWeight);
735  if ( histogram_tau2Pt && histogram_tau2Eta && histogram_tau2Phi ) {
736  histogram_tau2Pt->Fill(genLeg2->pt(), evtWeight);
737  histogram_tau2Eta->Fill(genLeg2->eta(), evtWeight);
738  histogram_tau2Phi->Fill(genLeg2->phi(), evtWeight);
739  }
740  if ( histogram_leg2Pt && histogram_leg2Eta && histogram_leg2Phi ) {
741  histogram_leg2Pt->Fill(visDecayProduct2->pt(), evtWeight);
742  histogram_leg2Eta->Fill(visDecayProduct2->eta(), evtWeight);
743  histogram_leg2Phi->Fill(visDecayProduct2->phi(), evtWeight);
744  }
745  histogram_leg2X->Fill(X2, evtWeight);
746  if ( histogram_leg2XforLeg1X0_00to0_25 && histogram_leg2XforLeg1X0_25to0_50 && histogram_leg2XforLeg1X0_50to0_75 && histogram_leg2XforLeg1X0_75to1_00 ) {
747  if ( X1 < 0.25 ) histogram_leg2XforLeg1X0_00to0_25->Fill(X2, evtWeight);
748  else if ( X1 < 0.50 ) histogram_leg2XforLeg1X0_25to0_50->Fill(X2, evtWeight);
749  else if ( X1 < 0.75 ) histogram_leg2XforLeg1X0_50to0_75->Fill(X2, evtWeight);
750  else histogram_leg2XforLeg1X0_75to1_00->Fill(X2, evtWeight);
751  if ( histogram_leg2Mt ) histogram_leg2Mt->Fill(compMt(visDecayProduct2->p4(), metP4), evtWeight);
752  }
753  histogram_VisDeltaPhiLeg1Leg2->Fill(normalizedPhi(visDecayProduct1->phi() - visDecayProduct2->phi()), evtWeight);
754  }
755  }
756  }
757  }
758  }
759 }
760 
761 namespace
762 {
763  std::string runLumiEventNumbers_to_string(const edm::Event& evt)
764  {
765  edm::RunNumber_t run_number = evt.id().run();
767  edm::EventNumber_t event_number = evt.id().event();
768  std::ostringstream retVal;
769  retVal << "Run = " << run_number << ", LS = " << ls_number << ", Event = " << event_number;
770  return retVal.str();
771  }
772 }
773 
775 {
776  if ( verbosity_ ) {
777  std::cout << "<MCEmbeddingValidationAnalyzer::analyze>:" << std::endl;
778  std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
779  }
780 
781  const reco::Candidate* replacedMuonPlus = 0;
782  const reco::Candidate* replacedMuonMinus = 0;
783  if ( srcReplacedMuons_.label() != "" ) {
784  std::vector<reco::CandidateBaseRef> replacedMuons = getSelMuons(evt, srcReplacedMuons_);
785  for ( std::vector<reco::CandidateBaseRef>::const_iterator replacedMuon = replacedMuons.begin();
786  replacedMuon != replacedMuons.end(); ++replacedMuon ) {
787  if ( (*replacedMuon)->charge() > +0.5 ) replacedMuonPlus = replacedMuon->get();
788  else if ( (*replacedMuon)->charge() < -0.5 ) replacedMuonMinus = replacedMuon->get();
789  }
790  }
791  if ( verbosity_ ) {
792  if ( replacedMuonPlus ) std::cout << "replacedMuonPlus: Pt = " << replacedMuonPlus->pt() << ", eta = " << replacedMuonPlus->eta() << ", phi = " << replacedMuonPlus->phi() << std::endl;
793  if ( replacedMuonMinus ) std::cout << "replacedMuonMinus: Pt = " << replacedMuonMinus->pt() << ", eta = " << replacedMuonMinus->eta() << ", phi = " << replacedMuonMinus->phi() << std::endl;
794  }
795 
797  edm::Handle<CandidateView> genDiTaus;
798  evt.getByLabel(srcGenDiTaus_, genDiTaus);
799 
800  const reco::Candidate* genTauPlus = 0;
801  const reco::Candidate* genTauMinus = 0;
802  for ( CandidateView::const_iterator genDiTau = genDiTaus->begin();
803  genDiTau != genDiTaus->end(); ++genDiTau ) {
804  const reco::CompositeCandidate* genDiTau_composite = dynamic_cast<const reco::CompositeCandidate*>(&(*genDiTau));
805  if ( !(genDiTau_composite && genDiTau_composite->numberOfDaughters() == 2) ) continue;
806  size_t numDaughters = genDiTau_composite->numberOfDaughters();
807  for ( size_t iDaughter = 0; iDaughter < numDaughters; ++iDaughter ) {
808  const reco::Candidate* daughter = genDiTau_composite->daughter(iDaughter);
809  if ( daughter->charge() > +0.5 ) genTauPlus = daughter;
810  else if ( daughter->charge() < -0.5 ) genTauMinus = daughter;
811  }
812  }
813  if ( verbosity_ ) {
814  if ( genTauPlus ) std::cout << "genTauPlus: Pt = " << genTauPlus->pt() << ", eta = " << genTauPlus->eta() << ", phi = " << genTauPlus->phi() << std::endl;
815  if ( genTauMinus ) std::cout << "genTauMinus: Pt = " << genTauMinus->pt() << ", eta = " << genTauMinus->eta() << ", phi = " << genTauMinus->phi() << std::endl;
816  }
817 
818  double dRmin = 1.e+3;
819  if ( replacedMuonPlus && genTauPlus ) {
820  double dR = deltaR(genTauPlus->p4(), replacedMuonPlus->p4());
821  if ( verbosity_ ) std::cout << " dR(replacedMuonPlus, genTauPlus) = " << dR << std::endl;
822  if ( dR < dRmin ) dRmin = dR;
823  }
824  if ( replacedMuonPlus && genTauMinus ) {
825  double dR = deltaR(genTauMinus->p4(), replacedMuonPlus->p4());
826  if ( verbosity_ ) std::cout << " dR(replacedMuonPlus, genTauMinus) = " << dR << std::endl;
827  if ( dR < dRmin ) dRmin = dR;
828  }
829  if ( replacedMuonMinus && genTauPlus ) {
830  double dR = deltaR(genTauPlus->p4(), replacedMuonMinus->p4());
831  if ( verbosity_ ) std::cout << " dR(replacedMuonMinus, genTauPlus) = " << dR << std::endl;
832  if ( dR < dRmin ) dRmin = dR;
833  }
834  if ( replacedMuonMinus && genTauMinus ) {
835  double dR = deltaR(genTauMinus->p4(), replacedMuonMinus->p4());
836  if ( verbosity_ ) std::cout << " dR(replacedMuonMinus, genTauMinus) = " << dR << std::endl;
837  if ( dR < dRmin ) dRmin = dR;
838  }
839  if ( verbosity_ ) std::cout << "--> dRmin = " << dRmin << std::endl;
840  if ( dRmin < dRminSeparation_ ) return;
841  if ( verbosity_ ) std::cout << " cut on dRminSeparation = " << dRminSeparation_ << " passed." << std::endl;
842 
843 //--- compute event weight
844  double evtWeight = 1.0;
845  std::map<std::string, double> evtWeightMap;
846  for ( vInputTag::const_iterator srcWeight = srcOtherWeights_.begin();
847  srcWeight != srcOtherWeights_.end(); ++srcWeight ) {
849  evt.getByLabel(*srcWeight, weight);
850  //std::cout << "weight(" << srcWeight->label().data() << ":" << srcWeight->instance().data() << ") = " << (*weight) << std::endl;
851  evtWeight *= (*weight);
852  evtWeightMap[Form("%s_%s", srcWeight->label().data(), srcWeight->instance().data())] = (*weight);
853  }
854  if ( srcGenFilterInfo_.label() != "" ) {
855  edm::Handle<GenFilterInfo> genFilterInfo;
856  evt.getByLabel(srcGenFilterInfo_, genFilterInfo);
857  //std::cout << "genFilterInfo: numEventsTried = " << genFilterInfo->numEventsTried() << ", numEventsPassed = " << genFilterInfo->numEventsPassed() << std::endl;
858  if ( genFilterInfo->numEventsTried() > 0 ) {
859  double weight = genFilterInfo->filterEfficiency();
860  //std::cout << "weight(genFilterInfo) = " << weight << std::endl;
861  histogramGenFilterEfficiency_->Fill(weight, evtWeight);
862  evtWeight *= weight;
863  evtWeightMap[Form("%s_%s", srcGenFilterInfo_.label().data(), srcGenFilterInfo_.instance().data())] = weight;
864  }
865  }
866 
867  if ( evtWeight < 1.e-3 || evtWeight > 1.e+3 || TMath::IsNaN(evtWeight) ) return;
868 
869  histogramEventCounter_->Fill(1., evtWeight);
870 
871  double muonRadCorrWeight = 1.;
872  double muonRadCorrWeightUp = 1.;
873  double muonRadCorrWeightDown = 1.;
875  edm::Handle<double> muonRadCorrWeight_handle;
876  evt.getByLabel(srcMuonRadCorrWeight_, muonRadCorrWeight_handle);
877  muonRadCorrWeight = (*muonRadCorrWeight_handle);
878  evtWeightMap["muonRadCorrWeight"] = muonRadCorrWeight;
879  edm::Handle<double> muonRadCorrWeightUp_handle;
880  evt.getByLabel(srcMuonRadCorrWeightUp_, muonRadCorrWeightUp_handle);
881  muonRadCorrWeightUp = (*muonRadCorrWeightUp_handle);
882  edm::Handle<double> muonRadCorrWeightDown_handle;
883  evt.getByLabel(srcMuonRadCorrWeightDown_, muonRadCorrWeightDown_handle);
884  muonRadCorrWeightDown = (*muonRadCorrWeightDown_handle);
885  }
886  //std::cout << " muonRadCorrWeight = " << muonRadCorrWeight
887  // << " + " << (muonRadCorrWeightUp - muonRadCorrWeight)
888  // << " - " << (muonRadCorrWeight - muonRadCorrWeightDown) << std::endl;
889 
890  for ( std::vector<plotEntryTypeEvtWeight*>::iterator evtWeightPlotEntry = evtWeightPlotEntries_.begin();
891  evtWeightPlotEntry != evtWeightPlotEntries_.end(); ++evtWeightPlotEntry ) {
892  (*evtWeightPlotEntry)->fillHistograms(evt, evtWeightMap);
893  }
894 
895 //--- fill all histograms
897  evt.getByLabel(srcRecTracks_, tracks);
898  int numTracksPtGt5 = 0;
899  int numTracksPtGt10 = 0;
900  int numTracksPtGt20 = 0;
901  int numTracksPtGt30 = 0;
902  int numTracksPtGt40 = 0;
903  for ( reco::TrackCollection::const_iterator track = tracks->begin();
904  track != tracks->end(); ++track ) {
905  if ( track->pt() > 5. ) ++numTracksPtGt5;
906  if ( track->pt() > 10. ) ++numTracksPtGt10;
907  if ( track->pt() > 20. ) ++numTracksPtGt20;
908  if ( track->pt() > 30. ) ++numTracksPtGt30;
909  if ( track->pt() > 40. ) ++numTracksPtGt40;
910  }
911  histogramNumTracksPtGt5_->Fill(numTracksPtGt5, muonRadCorrWeight*evtWeight);
912  histogramNumTracksPtGt10_->Fill(numTracksPtGt10, muonRadCorrWeight*evtWeight);
913  histogramNumTracksPtGt20_->Fill(numTracksPtGt20, muonRadCorrWeight*evtWeight);
914  histogramNumTracksPtGt30_->Fill(numTracksPtGt30, muonRadCorrWeight*evtWeight);
915  histogramNumTracksPtGt40_->Fill(numTracksPtGt40, muonRadCorrWeight*evtWeight);
916 
918  evt.getByLabel(srcRecPFCandidates_, pfCandidates);
919  int numChargedPFCandsPtGt5 = 0;
920  int numChargedPFCandsPtGt10 = 0;
921  int numChargedPFCandsPtGt20 = 0;
922  int numChargedPFCandsPtGt30 = 0;
923  int numChargedPFCandsPtGt40 = 0;
924  int numNeutralPFCandsPtGt5 = 0;
925  int numNeutralPFCandsPtGt10 = 0;
926  int numNeutralPFCandsPtGt20 = 0;
927  int numNeutralPFCandsPtGt30 = 0;
928  int numNeutralPFCandsPtGt40 = 0;
929  for ( reco::PFCandidateCollection::const_iterator pfCandidate = pfCandidates->begin();
930  pfCandidate != pfCandidates->end(); ++pfCandidate ) {
931  if ( TMath::Abs(pfCandidate->charge()) > 0.5 ) {
932  if ( pfCandidate->pt() > 5. ) ++numChargedPFCandsPtGt5;
933  if ( pfCandidate->pt() > 10. ) ++numChargedPFCandsPtGt10;
934  if ( pfCandidate->pt() > 20. ) ++numChargedPFCandsPtGt20;
935  if ( pfCandidate->pt() > 30. ) ++numChargedPFCandsPtGt30;
936  if ( pfCandidate->pt() > 40. ) ++numChargedPFCandsPtGt40;
937  } else {
938  if ( pfCandidate->pt() > 5. ) ++numNeutralPFCandsPtGt5;
939  if ( pfCandidate->pt() > 10. ) ++numNeutralPFCandsPtGt10;
940  if ( pfCandidate->pt() > 20. ) ++numNeutralPFCandsPtGt20;
941  if ( pfCandidate->pt() > 30. ) ++numNeutralPFCandsPtGt30;
942  if ( pfCandidate->pt() > 40. ) ++numNeutralPFCandsPtGt40;
943  }
944  }
945  histogramNumChargedPFCandsPtGt5_->Fill(numChargedPFCandsPtGt5, muonRadCorrWeight*evtWeight);
946  histogramNumChargedPFCandsPtGt10_->Fill(numChargedPFCandsPtGt10, muonRadCorrWeight*evtWeight);
947  histogramNumChargedPFCandsPtGt20_->Fill(numChargedPFCandsPtGt20, muonRadCorrWeight*evtWeight);
948  histogramNumChargedPFCandsPtGt30_->Fill(numChargedPFCandsPtGt30, muonRadCorrWeight*evtWeight);
949  histogramNumChargedPFCandsPtGt40_->Fill(numChargedPFCandsPtGt40, muonRadCorrWeight*evtWeight);
950  histogramNumNeutralPFCandsPtGt5_->Fill(numNeutralPFCandsPtGt5, muonRadCorrWeight*evtWeight);
951  histogramNumNeutralPFCandsPtGt10_->Fill(numNeutralPFCandsPtGt10, muonRadCorrWeight*evtWeight);
952  histogramNumNeutralPFCandsPtGt20_->Fill(numNeutralPFCandsPtGt20, muonRadCorrWeight*evtWeight);
953  histogramNumNeutralPFCandsPtGt30_->Fill(numNeutralPFCandsPtGt30, muonRadCorrWeight*evtWeight);
954  histogramNumNeutralPFCandsPtGt40_->Fill(numNeutralPFCandsPtGt40, muonRadCorrWeight*evtWeight);
955 
957  evt.getByLabel(srcRecJets_, jets);
958  int numJetsRawPtGt20 = 0;
959  int numJetsRawPtGt20AbsEtaLt2_5 = 0;
960  int numJetsRawPtGt20AbsEta2_5to4_5 = 0;
961  int numJetsCorrPtGt20 = 0;
962  int numJetsCorrPtGt20AbsEtaLt2_5 = 0;
963  int numJetsCorrPtGt20AbsEta2_5to4_5 = 0;
964  int numJetsRawPtGt30 = 0;
965  int numJetsRawPtGt30AbsEtaLt2_5 = 0;
966  int numJetsRawPtGt30AbsEta2_5to4_5 = 0;
967  int numJetsCorrPtGt30 = 0;
968  int numJetsCorrPtGt30AbsEtaLt2_5 = 0;
969  int numJetsCorrPtGt30AbsEta2_5to4_5 = 0;
970  for ( pat::JetCollection::const_iterator jet = jets->begin();
971  jet != jets->end(); ++jet ) {
972 
973  reco::Candidate::LorentzVector rawJetP4 = jet->correctedP4("Uncorrected");
974  double rawJetPt = rawJetP4.pt();
975  double rawJetAbsEta = TMath::Abs(rawJetP4.eta());
976  if ( rawJetAbsEta < 4.5 ) { // CV: do not consider any jet reconstructed outside eta range used in H -> tautau analysis
977  histogramRawJetPt_->Fill(rawJetPt, muonRadCorrWeight*evtWeight);
978  if ( rawJetAbsEta < 2.5 ) histogramRawJetPtAbsEtaLt2_5_->Fill(rawJetPt, muonRadCorrWeight*evtWeight);
979  else if ( rawJetAbsEta < 4.5 ) histogramRawJetPtAbsEta2_5to4_5_->Fill(rawJetPt, muonRadCorrWeight*evtWeight);
980  if ( rawJetPt > 20. ) {
981  histogramRawJetEtaPtGt20_->Fill(rawJetP4.eta(), muonRadCorrWeight*evtWeight);
982  ++numJetsRawPtGt20;
983  if ( rawJetAbsEta < 2.5 ) ++numJetsRawPtGt20AbsEtaLt2_5;
984  else if ( rawJetAbsEta < 4.5 ) ++numJetsRawPtGt20AbsEta2_5to4_5;
985  }
986  if ( rawJetPt > 30. ) {
987  histogramRawJetEtaPtGt20_->Fill(rawJetP4.eta(), muonRadCorrWeight*evtWeight);
988  ++numJetsRawPtGt30;
989  if ( rawJetAbsEta < 2.5 ) ++numJetsRawPtGt30AbsEtaLt2_5;
990  else if ( rawJetAbsEta < 4.5 ) ++numJetsRawPtGt30AbsEta2_5to4_5;
991  }
992  }
993 
994  reco::Candidate::LorentzVector corrJetP4 = jet->p4();
995  double corrJetPt = corrJetP4.pt();
996  double corrJetAbsEta = TMath::Abs(corrJetP4.eta());
997  if ( corrJetAbsEta < 4.5 ) { // CV: do not consider any jet reconstructed outside eta range used in H -> tautau analysis
998  histogramCorrJetPt_->Fill(corrJetPt, muonRadCorrWeight*evtWeight);
999  if ( corrJetAbsEta < 2.5 ) histogramCorrJetPtAbsEtaLt2_5_->Fill(corrJetPt, muonRadCorrWeight*evtWeight);
1000  else if ( corrJetAbsEta < 4.5 ) histogramCorrJetPtAbsEta2_5to4_5_->Fill(corrJetPt, muonRadCorrWeight*evtWeight);
1001  if ( corrJetPt > 20. ) {
1002  histogramCorrJetEtaPtGt20_->Fill(corrJetP4.eta(), muonRadCorrWeight*evtWeight);
1003  ++numJetsCorrPtGt20;
1004  if ( corrJetAbsEta < 2.5 ) ++numJetsCorrPtGt20AbsEtaLt2_5;
1005  else if ( corrJetAbsEta < 4.5 ) ++numJetsCorrPtGt20AbsEta2_5to4_5;
1006  }
1007  if ( corrJetPt > 30. ) {
1008  histogramCorrJetEtaPtGt20_->Fill(corrJetP4.eta(), muonRadCorrWeight*evtWeight);
1009  ++numJetsCorrPtGt30;
1010  if ( corrJetAbsEta < 2.5 ) ++numJetsCorrPtGt30AbsEtaLt2_5;
1011  else if ( corrJetAbsEta < 4.5 ) ++numJetsCorrPtGt30AbsEta2_5to4_5;
1012  }
1013  }
1014  }
1015  histogramNumJetsRawPtGt20_->Fill(numJetsRawPtGt20, muonRadCorrWeight*evtWeight);
1016  histogramNumJetsRawPtGt20AbsEtaLt2_5_->Fill(numJetsRawPtGt20AbsEtaLt2_5, muonRadCorrWeight*evtWeight);
1017  histogramNumJetsRawPtGt20AbsEta2_5to4_5_->Fill(numJetsRawPtGt20AbsEta2_5to4_5, muonRadCorrWeight*evtWeight);
1018  histogramNumJetsCorrPtGt20_->Fill(numJetsCorrPtGt20, muonRadCorrWeight*evtWeight);
1019  histogramNumJetsCorrPtGt20AbsEtaLt2_5_->Fill(numJetsCorrPtGt20AbsEtaLt2_5, muonRadCorrWeight*evtWeight);
1020  histogramNumJetsCorrPtGt20AbsEta2_5to4_5_->Fill(numJetsCorrPtGt20AbsEta2_5to4_5, muonRadCorrWeight*evtWeight);
1021  histogramNumJetsRawPtGt30_->Fill(numJetsRawPtGt30, muonRadCorrWeight*evtWeight);
1022  histogramNumJetsRawPtGt30AbsEtaLt2_5_->Fill(numJetsRawPtGt30AbsEtaLt2_5, muonRadCorrWeight*evtWeight);
1023  histogramNumJetsRawPtGt30AbsEta2_5to4_5_->Fill(numJetsRawPtGt30AbsEta2_5to4_5, muonRadCorrWeight*evtWeight);
1024  histogramNumJetsCorrPtGt30_->Fill(numJetsCorrPtGt30, muonRadCorrWeight*evtWeight);
1025  histogramNumJetsCorrPtGt30AbsEtaLt2_5_->Fill(numJetsCorrPtGt30AbsEtaLt2_5, muonRadCorrWeight*evtWeight);
1026  histogramNumJetsCorrPtGt30AbsEta2_5to4_5_->Fill(numJetsCorrPtGt30AbsEta2_5to4_5, muonRadCorrWeight*evtWeight);
1027 
1029  evt.getByLabel(srcRecMuons_, muons);
1030  int numGlobalMuons = 0;
1031  int numStandAloneMuons = 0;
1032  int numPFMuons = 0;
1033  for ( reco::MuonCollection::const_iterator muon = muons->begin();
1034  muon != muons->end(); ++muon ) {
1035  if ( muon->isGlobalMuon() ) ++numGlobalMuons;
1036  if ( muon->isStandAloneMuon() ) ++numStandAloneMuons;
1037  if ( muon->isPFMuon() ) ++numPFMuons;
1038  }
1039  histogramNumGlobalMuons_->Fill(numGlobalMuons, muonRadCorrWeight*evtWeight);
1040  histogramNumStandAloneMuons_->Fill(numStandAloneMuons, muonRadCorrWeight*evtWeight);
1041  histogramNumPFMuons_->Fill(numPFMuons, muonRadCorrWeight*evtWeight);
1042 
1044  evt.getByLabel(srcTheRecVertex_, theVertex);
1045  if ( theVertex->size() >= 1 ) {
1046  const reco::Vertex::Point& theVertexPosition = theVertex->front().position();
1047  histogramTheRecVertexX_->Fill(theVertexPosition.x(), muonRadCorrWeight*evtWeight);
1048  histogramTheRecVertexY_->Fill(theVertexPosition.y(), muonRadCorrWeight*evtWeight);
1049  histogramTheRecVertexZ_->Fill(theVertexPosition.z(), muonRadCorrWeight*evtWeight);
1050  }
1052  evt.getByLabel(srcRecVertices_, vertices);
1053  for ( reco::VertexCollection::const_iterator vertex = vertices->begin();
1054  vertex != vertices->end(); ++vertex ) {
1055  histogramRecVertexX_->Fill(vertex->position().x(), muonRadCorrWeight*evtWeight);
1056  histogramRecVertexY_->Fill(vertex->position().y(), muonRadCorrWeight*evtWeight);
1057  histogramRecVertexZ_->Fill(vertex->position().z(), muonRadCorrWeight*evtWeight);
1058  }
1059  histogramNumRecVertices_->Fill(vertices->size(), muonRadCorrWeight*evtWeight);
1060  edm::Handle<reco::VertexCollection> verticesWithBS;
1061  evt.getByLabel(srcRecVerticesWithBS_, verticesWithBS);
1062  for ( reco::VertexCollection::const_iterator vertex = verticesWithBS->begin();
1063  vertex != verticesWithBS->end(); ++vertex ) {
1064  histogramRecVertexWithBSx_->Fill(vertex->position().x(), muonRadCorrWeight*evtWeight);
1065  histogramRecVertexWithBSy_->Fill(vertex->position().y(), muonRadCorrWeight*evtWeight);
1066  histogramRecVertexWithBSz_->Fill(vertex->position().z(), muonRadCorrWeight*evtWeight);
1067  }
1068  histogramNumRecVerticesWithBS_->Fill(verticesWithBS->size(), muonRadCorrWeight*evtWeight);
1069 
1071  evt.getByLabel(srcBeamSpot_, beamSpot);
1072  if ( beamSpot.isValid() ) {
1073  histogramBeamSpotX_->Fill(beamSpot->position().x(), muonRadCorrWeight*evtWeight);
1074  histogramBeamSpotY_->Fill(beamSpot->position().y(), muonRadCorrWeight*evtWeight);
1075  }
1076 
1077  for ( CandidateView::const_iterator genDiTau = genDiTaus->begin();
1078  genDiTau != genDiTaus->end(); ++genDiTau ) {
1079  histogramGenDiTauPt_->Fill(genDiTau->pt(), muonRadCorrWeight*evtWeight);
1080  histogramGenDiTauEta_->Fill(genDiTau->eta(), muonRadCorrWeight*evtWeight);
1081  histogramGenDiTauPhi_->Fill(genDiTau->phi(), muonRadCorrWeight*evtWeight);
1082  histogramGenDiTauMass_->Fill(genDiTau->mass(), muonRadCorrWeight*evtWeight);
1083  }
1084 
1085  bool passesCutsBeforeRotation = false;
1086  if ( (replacedMuonPlus && replacedMuonPlus->pt() > replacedMuonPtThresholdHigh_ && replacedMuonMinus && replacedMuonMinus->pt() > replacedMuonPtThresholdLow_) ||
1087  (replacedMuonMinus && replacedMuonMinus->pt() > replacedMuonPtThresholdHigh_ && replacedMuonPlus && replacedMuonPlus->pt() > replacedMuonPtThresholdLow_) ) passesCutsBeforeRotation = true;
1088  bool passesCutsAfterRotation = false;
1089  for ( CandidateView::const_iterator genDiTau = genDiTaus->begin();
1090  genDiTau != genDiTaus->end(); ++genDiTau ) {
1091  const reco::CompositeCandidate* genDiTau_composite = dynamic_cast<const reco::CompositeCandidate*>(&(*genDiTau));
1092  if ( !(genDiTau_composite && genDiTau_composite->numberOfDaughters() == 2) ) continue;
1093  const reco::Candidate* genTau1 = genDiTau_composite->daughter(0);
1094  const reco::Candidate* genTau2 = genDiTau_composite->daughter(1);
1095  if ( !(genTau1 && genTau2) ) continue;
1096  if ( (genTau1->pt() > replacedMuonPtThresholdHigh_ && genTau2->pt() > replacedMuonPtThresholdLow_ ) ||
1097  (genTau1->pt() > replacedMuonPtThresholdLow_ && genTau2->pt() > replacedMuonPtThresholdHigh_) ) {
1098  passesCutsAfterRotation = true;
1099  break;
1100  }
1101  histogramRotationAngleMatrix_->Fill(passesCutsBeforeRotation, passesCutsAfterRotation, muonRadCorrWeight*evtWeight);
1102  }
1103  if ( genTauPlus && genTauMinus ) {
1104  histogramGenDeltaPhiLeg1Leg2_->Fill(normalizedPhi(genTauPlus->phi() - genTauMinus->phi()), muonRadCorrWeight*evtWeight);
1105  reco::Particle::LorentzVector diTauP4_lab = genTauPlus->p4() + genTauMinus->p4();
1106  ROOT::Math::Boost boost_to_rf(diTauP4_lab.BoostToCM());
1107  reco::Particle::LorentzVector diTauP4_rf = boost_to_rf(diTauP4_lab);
1108  reco::Particle::LorentzVector tauPlusP4_rf = boost_to_rf(genTauPlus->p4());
1109  if ( (diTauP4_rf.P()*tauPlusP4_rf.P()) > 0. ) {
1110  double cosGjAngle = (diTauP4_rf.px()*tauPlusP4_rf.px() + diTauP4_rf.py()*tauPlusP4_rf.py() + diTauP4_rf.pz()*tauPlusP4_rf.pz())/(diTauP4_rf.P()*tauPlusP4_rf.P());
1111  double gjAngle = TMath::ACos(cosGjAngle);
1112  histogramGenDiTauDecayAngle_->Fill(gjAngle, muonRadCorrWeight*evtWeight);
1113  }
1114  }
1115  if ( replacedMuonPlus && genTauPlus && replacedMuonMinus && genTauMinus ) {
1116  histogramRotationLegPlusDeltaR_->Fill(deltaR(genTauPlus->p4(), replacedMuonPlus->p4()), muonRadCorrWeight*evtWeight);
1117  histogramRotationLegMinusDeltaR_->Fill(deltaR(genTauMinus->p4(), replacedMuonMinus->p4()), muonRadCorrWeight*evtWeight);
1118 
1119  reco::Particle::LorentzVector diTauP4_lab = genTauPlus->p4() + genTauMinus->p4();
1120  histogramPhiRotLegPlus_->Fill(SVfit_namespace::phiLabFromLabMomenta(replacedMuonPlus->p4(), diTauP4_lab), muonRadCorrWeight*evtWeight);
1121  histogramPhiRotLegMinus_->Fill(SVfit_namespace::phiLabFromLabMomenta(replacedMuonMinus->p4(), diTauP4_lab), muonRadCorrWeight*evtWeight);
1122  }
1123 
1124  fillVisPtEtaPhiMassDistributions(evt, srcGenLeg1_, srcGenLeg2_, histogramGenVisDiTauPt_, histogramGenVisDiTauEta_, histogramGenVisDiTauPhi_, histogramGenVisDiTauMass_, muonRadCorrWeight*evtWeight);
1125  fillVisPtEtaPhiMassDistributions(evt, srcRecLeg1_, srcRecLeg2_, histogramRecVisDiTauPt_, histogramRecVisDiTauEta_, histogramRecVisDiTauPhi_, histogramRecVisDiTauMass_, muonRadCorrWeight*evtWeight);
1126 
1127  typedef edm::View<reco::MET> METView;
1128  edm::Handle<METView> genCaloMETs;
1129  evt.getByLabel(srcGenCaloMEt_, genCaloMETs);
1130  const reco::Candidate::LorentzVector& genCaloMEtP4 = genCaloMETs->front().p4();
1131  histogramGenCaloMEt_->Fill(genCaloMEtP4.pt(), muonRadCorrWeight*evtWeight);
1132  edm::Handle<METView> genPFMETs;
1133  evt.getByLabel(srcGenPFMEt_, genPFMETs);
1134  const reco::Candidate::LorentzVector& genPFMEtP4 = genPFMETs->front().p4();
1135  histogramGenPFMEt_->Fill(genPFMEtP4.pt(), muonRadCorrWeight*evtWeight);
1136 
1137  fillX1andX2Distributions(evt, srcGenDiTaus_, srcGenLeg1_, srcGenLeg2_, genPFMEtP4,
1146  histogramGenVisDeltaPhiLeg1Leg2_, muonRadCorrWeight*evtWeight);
1147 
1148  edm::Handle<METView> recPFMETs;
1149  evt.getByLabel(srcRecPFMEt_, recPFMETs);
1150  const reco::Candidate::LorentzVector& recPFMEtP4 = recPFMETs->front().p4();
1151 
1152  fillX1andX2Distributions(evt, srcGenDiTaus_, srcRecLeg1_, srcRecLeg2_, recPFMEtP4,
1153  0, 0, 0,
1154  0, 0, 0, histogramRecLeg1X_,
1155  0, 0, 0, 0,
1157  0, 0, 0,
1158  0, 0, 0, histogramRecLeg2X_,
1159  0, 0, 0, 0, histogramRecLeg2PFMt_,
1160  histogramRecVisDeltaPhiLeg1Leg2_, muonRadCorrWeight*evtWeight);
1161 
1163  evt.getByLabel(srcGenParticles_, genParticles);
1164  reco::Candidate::LorentzVector sumGenParticleP4;
1165  reco::Candidate::LorentzVector sumGenParticleP4_charged;
1166  for ( reco::GenParticleCollection::const_iterator genParticle = genParticles->begin();
1167  genParticle != genParticles->end(); ++genParticle ) {
1168  if ( genParticle->status() == 1 ) {
1169  int absPdgId = TMath::Abs(genParticle->pdgId());
1170  if ( absPdgId == 12 || absPdgId == 14 || absPdgId == 16 ) continue;
1171  sumGenParticleP4 += genParticle->p4();
1172  if ( TMath::Abs(genParticle->charge()) > 0.5 ) sumGenParticleP4_charged += genParticle->p4();
1173  }
1174  }
1175  histogramSumGenParticlePt_->Fill(sumGenParticleP4.pt(), muonRadCorrWeight*evtWeight);
1176 
1178  evt.getByLabel(srcCaloTowers_, caloTowers);
1179  reco::Candidate::LorentzVector sumCaloTowerP4_ecal;
1180  double sumEtCaloTowersECAL = 0.;
1181  reco::Candidate::LorentzVector sumCaloTowerP4_hcal;
1182  double sumEtCaloTowersHCAL = 0.;
1183  reco::Candidate::LorentzVector sumCaloTowerP4_hf;
1184  double sumEtCaloTowersHF = 0.;
1185  reco::Candidate::LorentzVector sumCaloTowerP4_ho;
1186  double sumEtCaloTowersHO = 0.;
1187  for ( CaloTowerCollection::const_iterator caloTower = caloTowers->begin();
1188  caloTower != caloTowers->end(); ++caloTower ) {
1189  if ( caloTower->energy() != 0. ) {
1190  double emFrac_ecal = caloTower->emEnergy()/caloTower->energy();
1191  double emFrac_hcal = (caloTower->energyInHB() + caloTower->energyInHE())/caloTower->energy();
1192  double emFrac_hf = caloTower->energyInHF()/caloTower->energy();
1193  double emFrac_ho = caloTower->energyInHO()/caloTower->energy();
1194  sumCaloTowerP4_ecal += (emFrac_ecal*caloTower->p4());
1195  sumEtCaloTowersECAL += (emFrac_ecal*caloTower->et());
1196  sumCaloTowerP4_hcal += (emFrac_hcal*caloTower->p4());
1197  sumEtCaloTowersHCAL += (emFrac_hcal*caloTower->et());
1198  sumCaloTowerP4_hf += (emFrac_hf*caloTower->p4());
1199  sumEtCaloTowersHF += (emFrac_hf*caloTower->et());
1200  sumCaloTowerP4_ho += (emFrac_ho*caloTower->p4());
1201  sumEtCaloTowersHO += (emFrac_ho*caloTower->et());
1202  }
1203  }
1204  histogramRecCaloMEtECAL_->Fill(sumCaloTowerP4_ecal.pt(), muonRadCorrWeight*evtWeight);
1205  histogramRecCaloSumEtECAL_->Fill(sumEtCaloTowersECAL, muonRadCorrWeight*evtWeight);
1206  histogramRecCaloMEtHCAL_->Fill(sumCaloTowerP4_hcal.pt(), muonRadCorrWeight*evtWeight);
1207  histogramRecCaloSumEtHCAL_->Fill(sumEtCaloTowersHCAL, muonRadCorrWeight*evtWeight);
1208  histogramRecCaloMEtHF_->Fill(sumCaloTowerP4_hf.pt(), muonRadCorrWeight*evtWeight);
1209  histogramRecCaloSumEtHF_->Fill(sumEtCaloTowersHF, muonRadCorrWeight*evtWeight);
1210  histogramRecCaloMEtHO_->Fill(sumCaloTowerP4_ho.pt(), muonRadCorrWeight*evtWeight);
1211  histogramRecCaloSumEtHO_->Fill(sumEtCaloTowersHO, muonRadCorrWeight*evtWeight);
1212 
1213  if ( srcMuonsBeforeRad_.label() != "" && srcMuonsAfterRad_.label() != "" ) {
1214  reco::Candidate::LorentzVector genMuonPlusP4_beforeRad;
1215  bool genMuonPlus_beforeRad_found = false;
1216  reco::Candidate::LorentzVector genMuonMinusP4_beforeRad;
1217  bool genMuonMinus_beforeRad_found = false;
1218  reco::Candidate::LorentzVector genMuonPlusP4_afterRad;
1219  bool genMuonPlus_afterRad_found = false;
1220  reco::Candidate::LorentzVector genMuonMinusP4_afterRad;
1221  bool genMuonMinus_afterRad_found = false;
1222 
1223  findMuons(evt, srcMuonsBeforeRad_, genMuonPlusP4_beforeRad, genMuonPlus_beforeRad_found, genMuonMinusP4_beforeRad, genMuonMinus_beforeRad_found);
1224  findMuons(evt, srcMuonsAfterRad_, genMuonPlusP4_afterRad, genMuonPlus_afterRad_found, genMuonMinusP4_afterRad, genMuonMinus_afterRad_found);
1225 
1226  bool genMuonPlus_found = (genMuonPlus_beforeRad_found && genMuonPlus_afterRad_found);
1227  bool genMuonMinus_found = (genMuonMinus_beforeRad_found && genMuonMinus_afterRad_found);
1228 
1229  if ( genTauPlus && genMuonPlus_found && genTauMinus && genMuonMinus_found ) {
1230  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator muonRadCorrUncertaintyPlotEntry = muonRadCorrUncertaintyPlotEntries_beforeRad_.begin();
1231  muonRadCorrUncertaintyPlotEntry != muonRadCorrUncertaintyPlotEntries_beforeRad_.end(); ++muonRadCorrUncertaintyPlotEntry ) {
1232  (*muonRadCorrUncertaintyPlotEntry)->fillHistograms(numJetsCorrPtGt30, genMuonPlusP4_beforeRad, genMuonMinusP4_beforeRad, evtWeight, muonRadCorrWeight, muonRadCorrWeightUp, muonRadCorrWeightDown);
1233  }
1234  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator muonRadCorrUncertaintyPlotEntry = muonRadCorrUncertaintyPlotEntries_afterRad_.begin();
1235  muonRadCorrUncertaintyPlotEntry != muonRadCorrUncertaintyPlotEntries_afterRad_.end(); ++muonRadCorrUncertaintyPlotEntry ) {
1236  (*muonRadCorrUncertaintyPlotEntry)->fillHistograms(numJetsCorrPtGt30, genMuonPlusP4_afterRad, genMuonMinusP4_afterRad, evtWeight, muonRadCorrWeight, muonRadCorrWeightUp, muonRadCorrWeightDown);
1237  }
1238  for ( std::vector<plotEntryTypeMuonRadCorrUncertainty*>::iterator muonRadCorrUncertaintyPlotEntry = muonRadCorrUncertaintyPlotEntries_afterRadAndCorr_.begin();
1239  muonRadCorrUncertaintyPlotEntry != muonRadCorrUncertaintyPlotEntries_afterRadAndCorr_.end(); ++muonRadCorrUncertaintyPlotEntry ) {
1240  (*muonRadCorrUncertaintyPlotEntry)->fillHistograms(numJetsCorrPtGt30, genTauPlus->p4(), genTauMinus->p4(), evtWeight, muonRadCorrWeight, muonRadCorrWeightUp, muonRadCorrWeightDown);
1241  }
1242  } else {
1244  edm::LogWarning ("<MCEmbeddingValidationAnalyzer::analyze>")
1245  << " (" << runLumiEventNumbers_to_string(evt) << ")" << std::endl
1246  << " Failed to match muons before and after radiation to embedded tau leptons !!" << std::endl;
1247  std::cout << "genTauPlus: ";
1248  if ( genTauPlus ) std::cout << "Pt = " << genTauPlus->pt() << ", eta = " << genTauPlus->eta() << ", phi = " << genTauPlus->phi() << std::endl;
1249  else std::cout << "NA" << std::endl;
1250  std::cout << "genMuonPlus (before Rad.): ";
1251  if ( genMuonPlus_found ) std::cout << "Pt = " << genMuonPlusP4_beforeRad.pt() << ", eta = " << genMuonPlusP4_beforeRad.eta() << ", phi = " << genMuonPlusP4_beforeRad.phi() << std::endl;
1252  else std::cout << "NA" << std::endl;
1253  std::cout << "genTauMinus: ";
1254  if ( genTauMinus ) std::cout << "Pt = " << genTauMinus->pt() << ", eta = " << genTauMinus->eta() << ", phi = " << genTauMinus->phi() << std::endl;
1255  else std::cout << "NA" << std::endl;
1256  std::cout << "genMuonMinus (before Rad.): ";
1257  if ( genMuonMinus_found ) std::cout << "Pt = " << genMuonMinusP4_beforeRad.pt() << ", eta = " << genMuonMinusP4_beforeRad.eta() << ", phi = " << genMuonMinusP4_beforeRad.phi() << std::endl;
1258  else std::cout << "NA" << std::endl;
1260  }
1261  }
1262  }
1263 
1265  evt.getByLabel(srcL1ETM_, l1METs);
1266  const reco::Candidate::LorentzVector& l1MEtP4 = l1METs->front().p4();
1267  //std::cout << "L1MEt: Pt = " << l1MEtP4.pt() << ", phi = " << l1MEtP4.phi() << " (Et = " << l1METs->front().etMiss() << ", Px = " << l1MEtP4.px() << ", Py = " << l1MEtP4.py() << ")" << std::endl;
1268  //if ( l1MEtP4.pt() > 75. ) std::cout << "--> CHECK !!" << std::endl;
1269  typedef edm::View<reco::MET> METView;
1270  edm::Handle<METView> recCaloMETs;
1271  evt.getByLabel(srcRecCaloMEt_, recCaloMETs);
1272  const reco::Candidate::LorentzVector& recCaloMEtP4 = recCaloMETs->front().p4();
1273  //std::cout << "recCaloMEt: Pt = " << recCaloMEtP4.pt() << ", phi = " << recCaloMEtP4.phi() << " (Px = " << recCaloMEtP4.px() << ", Py = " << recCaloMEtP4.py() << ")" << std::endl;
1274  for ( CandidateView::const_iterator genDiTau = genDiTaus->begin();
1275  genDiTau != genDiTaus->end(); ++genDiTau ) {
1276  const reco::CompositeCandidate* genDiTau_composite = dynamic_cast<const reco::CompositeCandidate*>(&(*genDiTau));
1277  if ( !(genDiTau_composite && genDiTau_composite->numberOfDaughters() == 2) ) continue;
1278  const reco::Candidate* genDaughter1 = genDiTau_composite->daughter(0);
1279  const reco::GenParticle* genTau1 = ( genDaughter1->hasMasterClone() ) ?
1280  dynamic_cast<const reco::GenParticle*>(&(*genDaughter1->masterClone())) : dynamic_cast<const reco::GenParticle*>(genDaughter1);
1281  const reco::Candidate* genDaughter2 = genDiTau_composite->daughter(1);
1282  const reco::GenParticle* genTau2 = ( genDaughter2->hasMasterClone() ) ?
1283  dynamic_cast<const reco::GenParticle*>(&(*genDaughter2->masterClone())) : dynamic_cast<const reco::GenParticle*>(genDaughter2);
1284  if ( !(genTau1 && genTau2) ) continue;
1285  std::string genTauDecayMode1 = getGenTauDecayMode(genTau1);
1286  std::string genTauDecayMode2 = getGenTauDecayMode(genTau2);
1287  std::string genTauDecayMode_ref;
1288  if ( genTauDecayMode1 == "electron" || genTauDecayMode1 == "muon" ) genTauDecayMode_ref = genTauDecayMode2;
1289  else if ( genTauDecayMode2 == "electron" || genTauDecayMode2 == "muon" ) genTauDecayMode_ref = genTauDecayMode1;
1290  for ( std::vector<plotEntryTypeL1ETM*>::iterator l1ETMplotEntry = l1ETMplotEntries_.begin();
1291  l1ETMplotEntry != l1ETMplotEntries_.end(); ++l1ETMplotEntry ) {
1292  (*l1ETMplotEntry)->fillHistograms(genTauDecayMode_ref, l1MEtP4, genCaloMEtP4, recCaloMEtP4, genDiTau->p4(), muonRadCorrWeight*evtWeight);
1293  }
1294  }
1295 
1296  fillHistograms(electronDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1297  fillHistograms(electronDistributionsExtra_, numJetsCorrPtGt30, evt, es, muonRadCorrWeight*evtWeight);
1298  fillHistograms(electronEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1299  fillHistograms(gsfElectronEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1300  fillHistograms(electronL1TriggerEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1301  fillHistograms(muonDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1302  fillHistograms(muonEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1303  fillHistograms(muonL1TriggerEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1304  fillHistograms(tauDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1305  fillHistograms(tauDistributionsExtra_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1306  fillHistograms(tauEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1307  fillHistograms(l1ElectronDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1308  fillHistograms(l1MuonDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1309  fillHistograms(l1TauDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1310  fillHistograms(l1CentralJetDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1311  fillHistograms(l1ForwardJetDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1312  fillHistograms(metDistributions_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1313  fillHistograms(metL1TriggerEfficiencies_, numJetsCorrPtGt30, evt, muonRadCorrWeight*evtWeight);
1314 
1315  // CV: Check for presence in the embedded event of high Pt tracks, charged PFCandidates and muons
1316  // reconstructed near the direction of the replaced muons
1317  // (indicating that maybe not all muon signals were removed in the embedding)
1318  std::vector<const reco::Candidate*> replacedMuons;
1319  if ( replacedMuonPlus ) replacedMuons.push_back(replacedMuonPlus);
1320  if ( replacedMuonMinus ) replacedMuons.push_back(replacedMuonMinus);
1321  for ( std::vector<const reco::Candidate*>::const_iterator replacedMuon = replacedMuons.begin();
1322  replacedMuon != replacedMuons.end(); ++replacedMuon ) {
1323  for ( reco::TrackCollection::const_iterator track = tracks->begin();
1324  track != tracks->end(); ++track ) {
1325  if ( track->pt() > 10. ) {
1326  double dR = deltaR(track->eta(), track->phi(), (*replacedMuon)->eta(), (*replacedMuon)->phi());
1327  if ( dR < 0.1 ) {
1328  edm::LogWarning("MCEmbeddingValidationAnalyzer")
1329  << " (" << runLumiEventNumbers_to_string(evt) << ")" << std::endl
1330  << " Found track: Pt = " << track->pt() << ", eta = " << track->eta() << ", phi = " << track->phi() << ", charge = " << track->charge() << " in direction of"
1331  << " a replaced muon: Pt = " << (*replacedMuon)->pt()<< ", eta = " << (*replacedMuon)->eta() << ", phi = " << (*replacedMuon)->phi() << ", charge = " << (*replacedMuon)->charge()
1332  << " (dR = " << dR << "). This may point to a problem in removing all muon signals in the Embedding." << std::endl;
1333  histogramWarning_recTrackNearReplacedMuon_->Fill(track->charge()*(*replacedMuon)->charge(), muonRadCorrWeight*evtWeight);
1334  }
1335  }
1336  }
1337  for ( reco::PFCandidateCollection::const_iterator pfCandidate = pfCandidates->begin();
1338  pfCandidate != pfCandidates->end(); ++pfCandidate ) {
1339  if ( pfCandidate->pt() > 10. && TMath::Abs(pfCandidate->charge()) > 0.5 ) {
1340  double dR = deltaR(pfCandidate->eta(), pfCandidate->phi(), (*replacedMuon)->eta(), (*replacedMuon)->phi());
1341  if ( dR < 0.1 ) {
1342  edm::LogWarning("MCEmbeddingValidationAnalyzer")
1343  << " (" << runLumiEventNumbers_to_string(evt) << ")" << std::endl
1344  << " Found charged PFCandidate: Pt = " << pfCandidate->pt() << ", eta = " << pfCandidate->eta() << ", phi = " << pfCandidate->phi() << ", charge = " << pfCandidate->charge() << " in direction of"
1345  << " a replaced muon: Pt = " << (*replacedMuon)->pt()<< ", eta = " << (*replacedMuon)->eta() << ", phi = " << (*replacedMuon)->phi() << ", charge = " << (*replacedMuon)->charge()
1346  << " (dR = " << dR << "). This may point to a problem in removing all muon signals in the Embedding." << std::endl;
1347  histogramWarning_recPFCandNearReplacedMuon_->Fill(pfCandidate->charge()*(*replacedMuon)->charge(), muonRadCorrWeight*evtWeight);
1348  }
1349  }
1350  }
1351  for ( reco::MuonCollection::const_iterator muon = muons->begin();
1352  muon != muons->end(); ++muon ) {
1353  if ( muon->pt() > 10. ) {
1354  double dR = deltaR(muon->eta(), muon->phi(), (*replacedMuon)->eta(), (*replacedMuon)->phi());
1355  if ( dR < 0.1 ) {
1356  edm::LogWarning("MCEmbeddingValidationAnalyzer")
1357  << " (" << runLumiEventNumbers_to_string(evt) << ")" << std::endl
1358  << " Found track: Pt = " << muon->pt() << ", eta = " << muon->eta() << ", phi = " << muon->phi() << ", charge = " << muon->charge() << " in direction of"
1359  << " a replaced muon: Pt = " << (*replacedMuon)->pt()<< ", eta = " << (*replacedMuon)->eta() << ", phi = " << (*replacedMuon)->phi() << ", charge = " << (*replacedMuon)->charge()
1360  << " (dR = " << dR << "). This may point to a problem in removing all muon signals in the Embedding." << std::endl;
1361  histogramWarning_recMuonNearReplacedMuon_->Fill(muon->charge()*(*replacedMuon)->charge(), muonRadCorrWeight*evtWeight);
1362  }
1363  }
1364  }
1365  }
1366 }
1367 
1369 
RunNumber_t run() const
Definition: EventID.h:42
tuple keyword
Definition: statics.py:7
const double Pi
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:44
std::string genTauDecayMode(const reco::CompositePtrCandidate &c)
Definition: JetMCTag.cc:81
std::vector< metL1TriggerEfficiencyType * > metL1TriggerEfficiencies_
void fillHistograms(std::vector< T * > collection, int numJets, const edm::Event &evt, double evtWeight)
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
void setupL1ExtraObjectDistribution(int, int, const edm::ParameterSet &, const std::string &, std::vector< l1ExtraObjectDistributionT< T > * > &)
virtual double energy() const =0
energy
std::vector< plotEntryTypeMuonRadCorrUncertainty * > muonRadCorrUncertaintyPlotEntries_afterRad_
std::vector< leptonEfficiencyT< pat::Electron > * > electronEfficiencies_
void analyze(const edm::Event &, const edm::EventSetup &)
math::XYZVector Vector
point in the space
Definition: Candidate.h:47
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:872
double phiLabFromLabMomenta(const reco::Candidate::LorentzVector &motherP4, const reco::Candidate::LorentzVector &visP4)
unsigned int EventNumber_t
Definition: EventID.h:30
virtual float eta() const =0
momentum pseudorapidity
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< ParameterSet > VParameterSet
Definition: ParameterSet.h:33
std::vector< leptonL1TriggerEfficiencyT1T2< pat::Electron, l1extra::L1EmParticle > * > electronL1TriggerEfficiencies_
MCEmbeddingValidationAnalyzer(const edm::ParameterSet &)
std::vector< CaloTower >::const_iterator const_iterator
bool exists(std::string const &parameterName) const
checks if a parameter exists
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
virtual float phi() const =0
momentum azimuthal angle
void setupMEtL1TriggerEfficiency(int, int, const edm::ParameterSet &, const std::string &, std::vector< metL1TriggerEfficiencyType * > &)
std::vector< plotEntryTypeL1ETM * > l1ETMplotEntries_
std::vector< tauDistributionExtra * > tauDistributionsExtra_
unsigned int LuminosityBlockNumber_t
Definition: EventID.h:31
virtual size_type numberOfDaughters() const
number of daughters
std::vector< plotEntryTypeMuonRadCorrUncertainty * > muonRadCorrUncertaintyPlotEntries_afterRadAndCorr_
void Fill(long long x)
std::vector< leptonDistributionT< pat::Electron > * > electronDistributions_
std::vector< PFCandidatePtr > pfCandidates(const PFJet &jet, int particleId, bool sort=true)
edm::View< reco::Candidate > CandidateView
virtual bool hasMasterClone() const =0
virtual float pt() const =0
transverse momentum
std::vector< edm::InputTag > vInputTag
reco::Candidate::Vector compCrossProduct(const reco::Candidate::Vector &p1, const reco::Candidate::Vector &p2)
void bookHistograms(std::vector< T * > collection, DQMStore &dqmStore)
void setupLeptonL1TriggerEfficiency(int, int, const edm::ParameterSet &, const std::string &, std::vector< leptonL1TriggerEfficiencyT1T2< T1, T2 > * > &)
std::vector< leptonEfficiencyT< reco::GsfElectron > * > gsfElectronEfficiencies_
vector< PseudoJet > jets
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
edm::View< reco::MET > METView
reco::Candidate::Vector normalize(const reco::Candidate::Vector &p)
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
virtual int charge() const =0
electric charge
std::vector< leptonEfficiencyT< pat::Muon > * > muonEfficiencies_
bool isValid() const
Definition: HandleBase.h:76
double p2[4]
Definition: TauolaWrapper.h:90
std::string dqmDirectory_full(const std::string &dqmSubDirectory)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
std::vector< plotEntryTypeEvtWeight * > evtWeightPlotEntries_
std::vector< l1ExtraObjectDistributionT< l1extra::L1JetParticle > * > l1CentralJetDistributions_
double compScalarProduct(const reco::Candidate::Vector &p1, const reco::Candidate::Vector &p2)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
std::vector< plotEntryTypeMuonRadCorrUncertainty * > muonRadCorrUncertaintyPlotEntries_beforeRad_
std::vector< l1ExtraObjectDistributionT< l1extra::L1JetParticle > * > l1TauDistributions_
tuple tracks
Definition: testEve_cfg.py:39
std::vector< metDistributionType * > metDistributions_
virtual const Candidate * daughter(size_type) const
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
void setupLeptonDistribution(int, int, const edm::ParameterSet &, const std::string &, std::vector< leptonDistributionT< T > * > &)
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:41
std::string const & label() const
Definition: InputTag.h:42
void cleanCollection(std::vector< T * > collection)
edm::EventID id() const
Definition: EventBase.h:56
double p1[4]
Definition: TauolaWrapper.h:89
void setupMEtDistribution(int, int, const edm::ParameterSet &, const std::string &, std::vector< metDistributionType * > &)
tuple muons
Definition: patZpeak.py:38
std::vector< leptonL1TriggerEfficiencyT1T2< pat::Muon, l1extra::L1MuonParticle > * > muonL1TriggerEfficiencies_
void setupTauDistributionExtra(int, int, const edm::ParameterSet &, const std::string &, std::vector< tauDistributionExtra * > &)
void setupElectronDistributionExtra(int, int, const edm::ParameterSet &, const std::string &, std::vector< electronDistributionExtra * > &)
T get() const
get a component
Definition: Candidate.h:219
std::vector< l1ExtraObjectDistributionT< l1extra::L1EmParticle > * > l1ElectronDistributions_
std::vector< l1ExtraObjectDistributionT< l1extra::L1MuonParticle > * > l1MuonDistributions_
tuple cout
Definition: gather_cfg.py:121
std::vector< leptonDistributionT< pat::Muon > * > muonDistributions_
unsigned int RunNumber_t
Definition: EventRange.h:32
void findMuons(const edm::Event &, const edm::InputTag &, reco::Candidate::LorentzVector &, bool &, reco::Candidate::LorentzVector &, bool &)
std::vector< reco::CandidateBaseRef > getSelMuons(const edm::Event &, const edm::InputTag &)
Definition: DDAxes.h:10
int weight
Definition: histoStyle.py:50
std::vector< electronDistributionExtra * > electronDistributionsExtra_
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:1000
double normalizedPhi(double phi)
Definition: normalizedPhi.cc:5
moduleLabel_(iConfig.getParameter< string >("@module_label"))
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:27
std::vector< leptonEfficiencyT< pat::Tau > * > tauEfficiencies_
std::string const & instance() const
Definition: InputTag.h:43
std::vector< l1ExtraObjectDistributionT< l1extra::L1JetParticle > * > l1ForwardJetDistributions_
void setupLeptonEfficiency(int, int, const edm::ParameterSet &, const std::string &, std::vector< leptonEfficiencyT< T > * > &)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:584
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
virtual const CandidateBaseRef & masterClone() const =0
std::vector< leptonDistributionT< pat::Tau > * > tauDistributions_