CMS 3D CMS Logo

ResolutionAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: ResolutionAnalyzer
4 // Class: ResolutionAnalyzer
5 //
13 //
14 // Original Author: Marco De Mattia
15 // Created: Thu Sep 11 12:16:00 CEST 2008
16 //
17 //
18 
19 // system include files
20 #include <memory>
21 #include <string>
22 #include <vector>
23 #include <CLHEP/Vector/LorentzVector.h>
24 
25 // user include files
37 #include "HepMC/GenEvent.h"
38 #include "HepMC/GenParticle.h"
39 #include "HepPDT/ParticleDataTable.hh"
40 #include "HepPDT/TableBuilder.hh"
41 #include "HepPDT/defs.h"
50 
51 #include "Histograms.h"
52 #include "MuScleFitUtils.h"
53 //
54 // class decleration
55 //
56 
58 public:
59  explicit ResolutionAnalyzer(const edm::ParameterSet&);
60  ~ResolutionAnalyzer() override;
61 
62 private:
63  void analyze(const edm::Event&, const edm::EventSetup&) override;
64  void endJob() override{};
65 
66  template <typename T>
67  std::vector<reco::LeafCandidate> fillMuonCollection(const std::vector<T>& tracks) {
68  std::vector<reco::LeafCandidate> muons;
69  typename std::vector<T>::const_iterator track;
70  for (track = tracks.begin(); track != tracks.end(); ++track) {
72  track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
74  if (debug_ > 0)
75  std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon
76  << ": initial value Pt = " << mu.Pt() << std::endl;
77  reco::LeafCandidate muon(track->charge(), mu);
78  // Store muon
79  // ----------
80  muons.push_back(muon);
81  }
82  return muons;
83  }
84 
86  void fillHistoMap();
88  void writeHistoMap();
91 
92  // ----------member data ---------------------------
93 
94  // Collections labels
95  // ------------------
97 
101  bool debug_;
102  std::map<std::string, Histograms*> mapHisto_;
103  TFile* outputFile_;
104 
108 
109  TString treeFileName_;
110  int32_t maxEvents_;
111 
112  double ptMax_;
113 
119 };
120 
121 //
122 // constructors and destructor
123 //
125  : theMuonLabel_(iConfig.getParameter<edm::InputTag>("MuonLabel")),
126  theMuonType_(iConfig.getParameter<int>("MuonType")),
127  theRootFileName_(iConfig.getUntrackedParameter<std::string>("OutputFileName")),
128  theCovariancesRootFileName_(iConfig.getUntrackedParameter<std::string>("InputFileName")),
129  debug_(iConfig.getUntrackedParameter<bool>("Debug")),
130  resonance_(iConfig.getParameter<bool>("Resonance")),
131  readCovariances_(iConfig.getParameter<bool>("ReadCovariances")),
132  treeFileName_(iConfig.getParameter<std::string>("InputTreeName")),
133  maxEvents_(iConfig.getParameter<int>("MaxEvents")),
134  ptMax_(iConfig.getParameter<double>("PtMax")) {
135  //now do what ever initialization is needed
136 
137  // Initial parameters values
138  // -------------------------
139  int resolFitType = iConfig.getParameter<int>("ResolFitType");
140  MuScleFitUtils::ResolFitType = resolFitType;
141  // MuScleFitUtils::resolutionFunction = resolutionFunctionArray[resolFitType];
143  // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType];
145 
146  MuScleFitUtils::parResol = iConfig.getParameter<std::vector<double> >("parResol");
147 
148  MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("ResFind");
149 
150  outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
151  outputFile_->cd();
152  fillHistoMap();
153 
154  eventCounter_ = 0;
155 }
156 
158  outputFile_->cd();
159  writeHistoMap();
160  outputFile_->Close();
161  std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
162 }
163 
164 //
165 // member functions
166 //
167 
168 // ------------ method called to for each event ------------
170  using namespace edm;
171 
172  std::cout << "starting" << std::endl;
173 
174  lorentzVector nullLorentzVector(0, 0, 0, 0);
175 
176  RootTreeHandler rootTreeHandler;
177  typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
178  MuonPairVector savedPairVector;
179  MuonPairVector genPairVector;
180 
181  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
182  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &evtRun, &genPairVector);
183  MuonPairVector::iterator savedPair = savedPairVector.begin();
184  MuonPairVector::iterator genPair = genPairVector.begin();
185  std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
186  for (; savedPair != savedPairVector.end(); ++savedPair, ++genPair) {
187  ++eventCounter_;
188 
189  if ((eventCounter_ % 10000) == 0) {
190  std::cout << "event = " << eventCounter_ << std::endl;
191  }
192 
193  lorentzVector recMu1(savedPair->first);
194  lorentzVector recMu2(savedPair->second);
195 
196  if (resonance_) {
197  // Histograms with genParticles characteristics
198  // --------------------------------------------
199 
200  reco::Particle::LorentzVector genMother(genPair->first + genPair->second);
201 
202  mapHisto_["GenMother"]->Fill(genMother);
203  mapHisto_["DeltaGenMotherMuons"]->Fill(genPair->first, genPair->second);
204  mapHisto_["GenMotherMuons"]->Fill(genPair->first);
205  mapHisto_["GenMotherMuons"]->Fill(genPair->second);
206 
207  // Match the reco muons with the gen and sim tracks
208  // ------------------------------------------------
209  if (checkDeltaR(genPair->first, recMu1)) {
210  mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Pt() + recMu1.Pt()) / genPair->first.Pt(), -1);
211  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Theta() + recMu1.Theta()), -1);
212  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(
213  recMu1,
214  (-cos(genPair->first.Theta()) / sin(genPair->first.Theta()) + cos(recMu1.Theta()) / sin(recMu1.Theta())),
215  -1);
216  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Eta() + recMu1.Eta()), -1);
217  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
218  mapHisto_["PhiResolutionGenVSMu"]->Fill(
219  recMu1, MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()), -1);
220  recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
221  deltaPtOverPt_->Fill((recMu1.Pt() - genPair->first.Pt()) / genPair->first.Pt());
222  if (fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2) {
223  recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
224  deltaPtOverPtForEta12_->Fill((recMu1.Pt() - genPair->first.Pt()) / genPair->first.Pt());
225  }
226  }
227  if (checkDeltaR(genPair->second, recMu2)) {
228  mapHisto_["PtResolutionGenVSMu"]->Fill(
229  recMu2, (-genPair->second.Pt() + recMu2.Pt()) / genPair->second.Pt(), +1);
230  mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2, (-genPair->second.Theta() + recMu2.Theta()), +1);
231  mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(
232  recMu2,
233  (-cos(genPair->second.Theta()) / sin(genPair->second.Theta()) + cos(recMu2.Theta()) / sin(recMu2.Theta())),
234  +1);
235  mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2, (-genPair->second.Eta() + recMu2.Eta()), +1);
236  // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
237  mapHisto_["PhiResolutionGenVSMu"]->Fill(
238  recMu2, MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()), +1);
239  recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
240  deltaPtOverPt_->Fill((recMu2.Pt() - genPair->second.Pt()) / genPair->second.Pt());
241  if (fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2) {
242  recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
243  deltaPtOverPtForEta12_->Fill((recMu2.Pt() - genPair->second.Pt()) / genPair->second.Pt());
244  }
245  }
246 
247  // Fill the mass resolution histograms
248  // -----------------------------------
249  // check if the recoMuons match the genMuons
250  // if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
251  if (genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
252  checkDeltaR(genPair->first, recMu1) && checkDeltaR(genPair->second, recMu2)) {
253  double recoMass = (recMu1 + recMu2).mass();
254  double genMass = (genPair->first + genPair->second).mass();
255  // first is always mu-, second is always mu+
256  mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
257 
258  // Fill the reconstructed resonance
259  reco::Particle::LorentzVector recoResonance(recMu1 + recMu2);
260  mapHisto_["RecoResonance"]->Fill(recoResonance);
261  mapHisto_["DeltaRecoResonanceMuons"]->Fill(recMu1, recMu2);
262  mapHisto_["RecoResonanceMuons"]->Fill(recMu1);
263  mapHisto_["RecoResonanceMuons"]->Fill(recMu2);
264 
265  // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
266  if (genMass != 0) {
267  // double diffMass = (recoMass - genMass)/genMass;
268  double diffMass = recoMass - genMass;
269  // Fill if for both muons
270  double pt1 = recMu1.pt();
271  double eta1 = recMu1.eta();
272  double pt2 = recMu2.pt();
273  double eta2 = recMu2.eta();
274  // This is to avoid nan
275  if (diffMass == diffMass) {
276  massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
277  massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
278  } else {
279  std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
280  }
281  // Fill with mass resolution from resolution function
282  double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
283  // The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
284  // We alos take the squared value, since var = sigma^2.
285  mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
286  mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
287  }
288 
289  // Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
290  mapHisto_["hFunctionResolPt"]->Fill(
291  recMu1,
292  MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
293  -1);
294  mapHisto_["hFunctionResolCotgTheta"]->Fill(
295  recMu1,
296  MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
297  -1);
298  mapHisto_["hFunctionResolPhi"]->Fill(
299  recMu1,
300  MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
301  -1);
302  mapHisto_["hFunctionResolPt"]->Fill(
303  recMu2,
304  MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
305  +1);
306  mapHisto_["hFunctionResolCotgTheta"]->Fill(
307  recMu2,
308  MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
309  +1);
310  mapHisto_["hFunctionResolPhi"]->Fill(
311  recMu2,
312  MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
313  +1);
314 
315  if (readCovariances_) {
316  // Compute mass error terms
317  // ------------------------
318  double mass = (recMu1 + recMu2).mass();
319  double pt1 = recMu1.Pt();
320  double phi1 = recMu1.Phi();
321  double eta1 = recMu1.Eta();
322  double theta1 = 2 * atan(exp(-eta1));
323  double pt2 = recMu2.Pt();
324  double phi2 = recMu2.Phi();
325  double eta2 = recMu2.Eta();
326  double theta2 = 2 * atan(exp(-eta2));
327  // Derivatives
328  double mMu2 = MuScleFitUtils::mMu2;
329  double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
330  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
331  pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
332  mass;
333  double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
334  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
335  pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
336  mass;
337  double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
338  double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
339  double dmdcotgth1 =
340  (pt1 * pt1 * cos(theta1) / sin(theta1) *
341  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
342  pt1 * pt2 * cos(theta2) / sin(theta2)) /
343  mass;
344  double dmdcotgth2 =
345  (pt2 * pt2 * cos(theta2) / sin(theta2) *
346  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
347  pt2 * pt1 * cos(theta1) / sin(theta1)) /
348  mass;
349 
350  // Multiplied by the pt here
351  // -------------------------
352  double dmdpt[2] = {dmdpt1 * recMu1.Pt(), dmdpt2 * recMu2.Pt()};
353  double dmdphi[2] = {dmdphi1, dmdphi2};
354  double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
355 
356  // Evaluate the single terms in the mass error expression
357 
358  reco::Particle::LorentzVector* recMu[2] = {&recMu1, &recMu2};
359  int charge[2] = {-1, +1};
360 
361  double fullMassRes = 0.;
362  double massRes1 = 0.;
363  double massRes2 = 0.;
364  double massRes3 = 0.;
365  double massRes4 = 0.;
366  double massRes5 = 0.;
367  double massRes6 = 0.;
368  double massResPtAndPt12 = 0.;
369 
370  for (int i = 0; i < 2; ++i) {
371  double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
372  double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
373  double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
374  double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
375  double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
376  double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
377 
378  double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
379  double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
380  double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
381  double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
382  double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
383  double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
384 
385  // ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
386  mapHisto_["MassResolutionPt"]->Fill(*(recMu[i]), ptVariance * std::pow(dmdpt[i], 2), charge[i]);
387  mapHisto_["MassResolutionCotgTheta"]->Fill(
388  *(recMu[i]), cotgThetaVariance * std::pow(dmdcotgth[i], 2), charge[i]);
389  mapHisto_["MassResolutionPhi"]->Fill(*(recMu[i]), phiVariance * std::pow(dmdphi[i], 2), charge[i]);
390  mapHisto_["MassResolutionPt-CotgTheta"]->Fill(
391  *(recMu[i]), 2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i], charge[i]);
392  mapHisto_["MassResolutionPt-Phi"]->Fill(*(recMu[i]), 2 * pt_phi * dmdpt[i] * dmdphi[i], charge[i]);
393  mapHisto_["MassResolutionCotgTheta-Phi"]->Fill(
394  *(recMu[i]), 2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i], charge[i]);
395 
396  mapHisto_["MassResolutionPt1-Pt2"]->Fill(*(recMu[i]), pt1_pt2 * dmdpt[0] * dmdpt[1], charge[i]);
397  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill(
398  *(recMu[i]), cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1], charge[i]);
399  mapHisto_["MassResolutionPhi1-Phi2"]->Fill(*(recMu[i]), phi1_phi2 * dmdphi[0] * dmdphi[1], charge[i]);
400  // This must be filled for both configurations: 12 and 21
401  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill(
402  *(recMu[i]), pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1], charge[i]);
403  mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill(
404  *(recMu[i]), pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0], charge[i]);
405  mapHisto_["MassResolutionPt12-Phi21"]->Fill(*(recMu[i]), pt12_phi21 * dmdpt[0] * dmdphi[1], charge[i]);
406  mapHisto_["MassResolutionPt12-Phi21"]->Fill(*(recMu[i]), pt12_phi21 * dmdpt[1] * dmdphi[0], charge[i]);
407  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill(
408  *(recMu[i]), cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1], charge[i]);
409  mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill(
410  *(recMu[i]), cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0], charge[i]);
411 
412  // Sigmas for comparison
413  mapHisto_["sigmaPtFromVariance"]->Fill(*(recMu[i]), sqrt(ptVariance), charge[i]);
414  mapHisto_["sigmaCotgThetaFromVariance"]->Fill(*(recMu[i]), sqrt(cotgThetaVariance), charge[i]);
415  mapHisto_["sigmaPhiFromVariance"]->Fill(*(recMu[i]), sqrt(phiVariance), charge[i]);
416 
417  // Pt term from function
418  mapHisto_["MassResolutionPtFromFunction"]->Fill(
419  *(recMu[i]),
421  (recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol)) *
422  std::pow(dmdpt[i], 2),
423  charge[i]);
424 
425  fullMassRes += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
426  phiVariance * std::pow(dmdphi[i], 2) +
427 
428  // These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
429  2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
430  2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
431 
432  pt1_pt2 * dmdpt[0] * dmdpt[1] + cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] +
433  phi1_phi2 * dmdphi[0] * dmdphi[1] +
434 
435  // These are filled twice, because of the two combinations
436  pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
437  pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
438  cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
439 
440  massRes1 += ptVariance * std::pow(dmdpt[i], 2);
441  massRes2 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2);
442  massRes3 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
443  phiVariance * std::pow(dmdphi[i], 2);
444  massRes4 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
445  phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
446  2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
447  2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i];
448  massRes5 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
449  phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
450  2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
451  2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
452  cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1];
453  massRes6 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
454  phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
455  2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
456  2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
457  cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1] +
458  pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
459  pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
460  cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
461 
462  massResPtAndPt12 += ptVariance * std::pow(dmdpt[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1];
463 
464  // Derivatives
465  mapHisto_["DerivativePt"]->Fill(*(recMu[i]), dmdpt[i], charge[i]);
466  mapHisto_["DerivativeCotgTheta"]->Fill(*(recMu[i]), dmdcotgth[i], charge[i]);
467  mapHisto_["DerivativePhi"]->Fill(*(recMu[i]), dmdphi[i], charge[i]);
468  }
469  // Fill the complete resolution function with covariance terms
470  mapHisto_["FullMassResolution"]->Fill(*(recMu[0]), fullMassRes, charge[0]);
471  mapHisto_["FullMassResolution"]->Fill(*(recMu[1]), fullMassRes, charge[1]);
472 
473  mapHisto_["MassRes1"]->Fill(*(recMu[0]), massRes1, charge[0]);
474  mapHisto_["MassRes1"]->Fill(*(recMu[1]), massRes1, charge[1]);
475  mapHisto_["MassRes2"]->Fill(*(recMu[0]), massRes2, charge[0]);
476  mapHisto_["MassRes2"]->Fill(*(recMu[1]), massRes2, charge[1]);
477  mapHisto_["MassRes3"]->Fill(*(recMu[0]), massRes3, charge[0]);
478  mapHisto_["MassRes3"]->Fill(*(recMu[1]), massRes3, charge[1]);
479  mapHisto_["MassRes4"]->Fill(*(recMu[0]), massRes4, charge[0]);
480  mapHisto_["MassRes4"]->Fill(*(recMu[1]), massRes4, charge[1]);
481  mapHisto_["MassRes5"]->Fill(*(recMu[0]), massRes5, charge[0]);
482  mapHisto_["MassRes5"]->Fill(*(recMu[1]), massRes5, charge[1]);
483  mapHisto_["MassRes6"]->Fill(*(recMu[0]), massRes6, charge[0]);
484  mapHisto_["MassRes6"]->Fill(*(recMu[1]), massRes6, charge[1]);
485  mapHisto_["MassResPtAndPt12"]->Fill(*(recMu[0]), massResPtAndPt12, charge[0]);
486  mapHisto_["MassResPtAndPt12"]->Fill(*(recMu[1]), massResPtAndPt12, charge[1]);
487  } else {
488  // Fill the covariances histograms
489  mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
490  }
491  }
492  } // end if resonance
493  }
494  // else {
495  //
496  // // Loop on the recMuons
497  // std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
498  // for ( ; recMuon!=muons.end(); ++recMuon ) {
499  // int charge = recMuon->charge();
500  //
501  // lorentzVector recMu(recMuon->p4());
502  //
503  // // Find the matching MC muon
504  // const HepMC::GenEvent* Evt = evtMC->GetEvent();
505  // //Loop on generated particles
506  // std::map<double, lorentzVector> genAssocMap;
507  // HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
508  // for( ; part!=Evt->particles_end(); ++part ) {
509  // if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
510  // lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
511  // (*part)->momentum().pz(),(*part)->momentum().e()));
512  //
513  // double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
514  // ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
515  //
516  // // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
517  // // Only in this case we consider it matching.
518  // if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
519  // }
520  // }
521  // // Take the closest in deltaR
522  // lorentzVector genMu(genAssocMap.begin()->second);
523  //
524  // // Histograms with genParticles characteristics
525  // // --------------------------------------------
526  //
527  // if(checkDeltaR(genMu,recMu)){
528  // mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
529  // mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
530  // mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
531  // +cos(recMu.Theta())/sin(recMu.Theta())),charge);
532  // mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
533  // mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
534  // }
535  // }
536 }
537 
539  outputFile_->cd();
540 
541  // Resonances
542  // If no Z is required, use a smaller mass range.
543  double minMass = 0.;
544  double maxMass = 200.;
545  if (MuScleFitUtils::resfind[0] != 1) {
546  maxMass = 30.;
547  }
548  mapHisto_["GenMother"] = new HParticle(outputFile_, "GenMother", minMass, maxMass);
549  mapHisto_["SimResonance"] = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
550  mapHisto_["RecoResonance"] = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
551 
552  // Resonance muons
553  mapHisto_["GenMotherMuons"] = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
554  mapHisto_["SimResonanceMuons"] = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
555  mapHisto_["RecoResonanceMuons"] = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
556 
557  // Deltas between resonance muons
558  mapHisto_["DeltaGenMotherMuons"] = new HDelta(outputFile_, "DeltaGenMotherMuons");
559  mapHisto_["DeltaSimResonanceMuons"] = new HDelta(outputFile_, "DeltaSimResonanceMuons");
560  mapHisto_["DeltaRecoResonanceMuons"] = new HDelta(outputFile_, "DeltaRecoResonanceMuons");
561 
562  // //Reconstructed muon kinematics
563  // //-----------------------------
564  // mapHisto_["hRecBestMu"] = new HParticle ("hRecBestMu");
565  // mapHisto_["hRecBestMu_Acc"] = new HParticle ("hRecBestMu_Acc");
566  // mapHisto_["hDeltaRecBestMu"] = new HDelta ("hDeltaRecBestMu");
567 
568  // mapHisto_["hRecBestRes"] = new HParticle ("hRecBestRes");
569  // mapHisto_["hRecBestRes_Acc"] = new HParticle ("hRecBestRes_Acc");
570  // mapHisto_["hRecBestResVSMu"] = new HMassVSPart ("hRecBestResVSMu");
571 
572  //Resolution VS muon kinematic
573  //----------------------------
574  mapHisto_["PtResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "PtResolutionGenVSMu");
575  mapHisto_["PtResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "PtResolutionSimVSMu");
576  mapHisto_["EtaResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "EtaResolutionGenVSMu");
577  mapHisto_["EtaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "EtaResolutionSimVSMu");
578  mapHisto_["ThetaResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "ThetaResolutionGenVSMu");
579  mapHisto_["ThetaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "ThetaResolutionSimVSMu");
580  mapHisto_["CotgThetaResolutionGenVSMu"] =
581  new HResolutionVSPart(outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
582  mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "CotgThetaResolutionSimVSMu");
583  mapHisto_["PhiResolutionGenVSMu"] =
584  new HResolutionVSPart(outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
585  mapHisto_["PhiResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "PhiResolutionSimVSMu");
586 
587  // Covariances between muons kinematic quantities
588  // ----------------------------------------------
589  double ptMax = ptMax_;
590 
591  // Mass resolution
592  // ---------------
593  mapHisto_["MassResolution"] = new HMassResolutionVSPart(outputFile_, "MassResolution");
594 
595  // mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
596 
597  // Mass resolution vs (pt, eta) of the muons from MC
598  massResolutionVsPtEta_ = new HCovarianceVSxy("Mass", "Mass", 100, 0., ptMax, 60, -3, 3);
599  // Mass resolution vs (pt, eta) of the muons from function
600  recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
601  recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
602  deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
603  deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
604 
605  // Muons resolutions from resolution functions
606  // -------------------------------------------
607  int totBinsY = 60;
608  mapHisto_["hFunctionResolMass"] = new HFunctionResolution(outputFile_, "hFunctionResolMass", ptMax, totBinsY);
609  mapHisto_["hFunctionResolPt"] = new HFunctionResolution(outputFile_, "hFunctionResolPt", ptMax, totBinsY);
610  mapHisto_["hFunctionResolCotgTheta"] =
611  new HFunctionResolution(outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
612  mapHisto_["hFunctionResolPhi"] = new HFunctionResolution(outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
613 
614  if (readCovariances_) {
615  // Covariances read from file. Used to compare the terms in the expression of mass error
616  mapHisto_["ReadCovariances"] = new HCovarianceVSParts(theCovariancesRootFileName_, "Covariance");
617 
618  // Variances
619  mapHisto_["MassResolutionPt"] = new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPt", ptMax);
620  mapHisto_["MassResolutionCotgTheta"] =
621  new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassCotgTheta", ptMax);
622  mapHisto_["MassResolutionPhi"] = new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPhi", ptMax);
623  // Covariances
624  mapHisto_["MassResolutionPt-CotgTheta"] =
625  new HFunctionResolution(outputFile_, "functionResolMassPt-CotgTheta", ptMax, totBinsY);
626  mapHisto_["MassResolutionPt-Phi"] =
627  new HFunctionResolution(outputFile_, "functionResolMassPt-Phi", ptMax, totBinsY);
628  mapHisto_["MassResolutionCotgTheta-Phi"] =
629  new HFunctionResolution(outputFile_, "functionResolMassCotgTheta-Phi", ptMax, totBinsY);
630  mapHisto_["MassResolutionPt1-Pt2"] =
631  new HFunctionResolution(outputFile_, "functionResolMassPt1-Pt2", ptMax, totBinsY);
632  mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] =
633  new HFunctionResolution(outputFile_, "functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
634  mapHisto_["MassResolutionPhi1-Phi2"] =
635  new HFunctionResolution(outputFile_, "functionResolMassPhi1-Phi2", ptMax, totBinsY);
636  mapHisto_["MassResolutionPt12-CotgTheta21"] =
637  new HFunctionResolution(outputFile_, "functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
638  mapHisto_["MassResolutionPt12-Phi21"] =
639  new HFunctionResolution(outputFile_, "functionResolMassPt12-Phi21", ptMax, totBinsY);
640  mapHisto_["MassResolutionCotgTheta12-Phi21"] =
641  new HFunctionResolution(outputFile_, "functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
642 
643  mapHisto_["sigmaPtFromVariance"] = new HFunctionResolution(outputFile_, "sigmaPtFromVariance", ptMax, totBinsY);
644  mapHisto_["sigmaCotgThetaFromVariance"] =
645  new HFunctionResolution(outputFile_, "sigmaCotgThetaFromVariance", ptMax, totBinsY);
646  mapHisto_["sigmaPhiFromVariance"] = new HFunctionResolution(outputFile_, "sigmaPhiFromVariance", ptMax, totBinsY);
647 
648  // Derivatives
649  mapHisto_["DerivativePt"] = new HFunctionResolution(outputFile_, "derivativePt", ptMax);
650  mapHisto_["DerivativeCotgTheta"] = new HFunctionResolution(outputFile_, "derivativeCotgTheta", ptMax);
651  mapHisto_["DerivativePhi"] = new HFunctionResolution(outputFile_, "derivativePhi", ptMax);
652 
653  // Pt term from function
654  mapHisto_["MassResolutionPtFromFunction"] =
655  new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPtFromFunction", ptMax);
656 
657  mapHisto_["FullMassResolution"] = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
658  mapHisto_["MassRes1"] = new HFunctionResolution(outputFile_, "massRes1", ptMax);
659  mapHisto_["MassRes2"] = new HFunctionResolution(outputFile_, "massRes2", ptMax);
660  mapHisto_["MassRes3"] = new HFunctionResolution(outputFile_, "massRes3", ptMax);
661  mapHisto_["MassRes4"] = new HFunctionResolution(outputFile_, "massRes4", ptMax);
662  mapHisto_["MassRes5"] = new HFunctionResolution(outputFile_, "massRes5", ptMax);
663  mapHisto_["MassRes6"] = new HFunctionResolution(outputFile_, "massRes6", ptMax);
664  mapHisto_["MassResPtAndPt12"] = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
665  } else {
666  mapHisto_["Covariances"] = new HCovarianceVSParts(outputFile_, "Covariance", ptMax);
667  }
668 }
669 
671  for (std::map<std::string, Histograms*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
672  histo++) {
673  (*histo).second->Write();
674  }
675  outputFile_->cd();
677  recoPtVsgenPt_->Write();
678  recoPtVsgenPtEta12_->Write();
679  deltaPtOverPt_->Write();
680  deltaPtOverPtForEta12_->Write();
681 }
682 
684  const reco::Particle::LorentzVector& recMu) {
685  //first is always mu-, second is always mu+
686  double deltaR =
687  sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
688  ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
689  if (deltaR < 0.01)
690  return true;
691  else if (debug_ > 0)
692  std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
693  << " DOES NOT MATCH with generated muon from resonance: " << std::endl
694  << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
695  return false;
696 }
697 
698 //define this as a plug-in
static double deltaPhiNoFabs(const double &phi1, const double &phi2)
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance c...
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
std::map< std::string, Histograms * > mapHisto_
std::string theRootFileName_
static std::vector< double > parResol
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void Fill(const double &x, const double &y, const double &a, const double &b) override
Definition: Histograms.h:1989
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, std::vector< std::pair< unsigned int, unsigned long long > > *evtRun, MuonPairVector *genPair=nullptr)
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:70
std::vector< reco::LeafCandidate > fillMuonCollection(const std::vector< T > &tracks)
HCovarianceVSxy * massResolutionVsPtEta_
void analyze(const edm::Event &, const edm::EventSetup &) override
void fillHistoMap()
Used to fill the map with the histograms needed.
int iEvent
Definition: GenABIO.cc:224
resolutionFunctionBase< std::vector< double > > * resolutionFunctionVecService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier when receiving a std::...
Definition: Functions.cc:92
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
void Write() override
Definition: Histograms.h:2020
T sqrt(T t)
Definition: SSEVec.h:19
void writeHistoMap()
Writes the histograms in the map.
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::string theCovariancesRootFileName_
bool checkDeltaR(const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recMu)
Returns true if the two particles have DeltaR < cut.
auto const & tracks
cannot be loose
void endJob() override
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
static int ResolFitType
static int goodmuon
static const double mMu2
HLT enums.
static double deltaPhi(const double &phi1, const double &phi2)
static std::vector< int > resfind
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
A set of histograms for resolution.
Definition: Histograms.h:1201
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
static resolutionFunctionBase< double * > * resolutionFunction
ResolutionAnalyzer(const edm::ParameterSet &)
edm::InputTag theMuonLabel_