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