CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
ErrorsPropagationAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: ErrorsPropagationAnalyzer
4 // Class: ErrorsPropagationAnalyzer
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 
24 #include "FWCore/Framework/interface/Frameworkfwd.h"
25 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
26 #include "FWCore/Framework/interface/Event.h"
27 #include "FWCore/Framework/interface/MakerMacros.h"
28 #include "FWCore/ParameterSet/interface/ParameterSet.h"
29 
30 #include <TH1D.h>
31 #include <TProfile.h>
32 #include <TString.h>
33 #include <TCanvas.h>
34 #include <TGraphAsymmErrors.h>
35 #include <TROOT.h>
36 
37 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
38 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
39 #include "MuScleFitUtils.h"
40 #include "MuonAnalysis/MomentumScaleCalibration/interface/SigmaPtDiff.h"
41 
42 //
43 // class declaration
44 //
45 
46 class ErrorsPropagationAnalyzer : public edm::one::EDAnalyzer<> {
47 public:
48  explicit ErrorsPropagationAnalyzer(const edm::ParameterSet&);
49  ~ErrorsPropagationAnalyzer() override;
50 
51 private:
52  void analyze(const edm::Event&, const edm::EventSetup&) override;
53  void fillHistograms();
54  void drawHistograms(const TProfile* histo,
55  const TProfile* histoPlusErr,
56  const TProfile* histoMinusErr,
57  const TString& type,
58  const TString& yLabel);
59  void fillValueError();
60  void endJob() override{};
62  double massResolution(const lorentzVector& mu1,
63  const lorentzVector& mu2,
64  const std::vector<double>& parval,
65  const double& sigmaPt1,
66  const double& sigmaPt2);
67  double massResolution(const lorentzVector& mu1,
68  const lorentzVector& mu2,
69  double* parval,
70  const double& sigmaPt1,
71  const double& sigmaPt2);
72 
73  TString treeFileName_;
74  int resolFitType_;
75  uint32_t maxEvents_;
76  TString outputFileName_;
77  int ptBins_;
78  double ptMin_;
79  double ptMax_;
80  int etaBins_;
81  double etaMin_;
82  double etaMax_;
83  bool debug_;
84 
85  double ptMinCut_, ptMaxCut_, etaMinCut_, etaMaxCut_;
86 
87  std::vector<double> parameters_;
88  std::vector<double> errors_;
89  std::vector<int> errorFactors_;
90 
91  std::vector<double> valuePlusError_;
92  std::vector<double> valueMinusError_;
93 
94  TProfile* sigmaPtVsEta_;
95  TProfile* sigmaPtVsEtaPlusErr_;
96  TProfile* sigmaPtVsEtaMinusErr_;
97 
98  TProfile* sigmaPtVsPt_;
99  TProfile* sigmaPtVsPtPlusErr_;
100  TProfile* sigmaPtVsPtMinusErr_;
101 
102  TProfile* sigmaPtVsEtaDiff_;
103  TProfile* sigmaPtVsPtDiff_;
104 
105  // Mass resolution
106  TProfile* sigmaMassVsEta_;
107  TProfile* sigmaMassVsEtaPlusErr_;
108  TProfile* sigmaMassVsEtaMinusErr_;
109 
110  TProfile* sigmaMassVsPt_;
111  TProfile* sigmaMassVsPtPlusErr_;
112  TProfile* sigmaMassVsPtMinusErr_;
113 
114  TProfile* sigmaMassOverMassVsEta_;
115  TProfile* sigmaMassOverMassVsEtaPlusErr_;
116  TProfile* sigmaMassOverMassVsEtaMinusErr_;
117 
118  TProfile* sigmaMassOverMassVsPt_;
119  TProfile* sigmaMassOverMassVsPtPlusErr_;
120  TProfile* sigmaMassOverMassVsPtMinusErr_;
121 };
122 
123 ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig)
124  : treeFileName_(iConfig.getParameter<std::string>("InputFileName")),
125  resolFitType_(iConfig.getParameter<int>("ResolFitType")),
126  maxEvents_(iConfig.getParameter<int>("MaxEvents")),
127  outputFileName_(iConfig.getParameter<std::string>("OutputFileName")),
128  ptBins_(iConfig.getParameter<int>("PtBins")),
129  ptMin_(iConfig.getParameter<double>("PtMin")),
130  ptMax_(iConfig.getParameter<double>("PtMax")),
131  etaBins_(iConfig.getParameter<int>("EtaBins")),
132  etaMin_(iConfig.getParameter<double>("EtaMin")),
133  etaMax_(iConfig.getParameter<double>("EtaMax")),
134  debug_(iConfig.getParameter<bool>("Debug")),
135  ptMinCut_(iConfig.getUntrackedParameter<double>("PtMinCut", 0.)),
136  ptMaxCut_(iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.)),
137  etaMinCut_(iConfig.getUntrackedParameter<double>("EtaMinCut", 0.)),
138  etaMaxCut_(iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.)) {
139  parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
140  errors_ = iConfig.getParameter<std::vector<double> >("Errors");
141  errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
142 
143  if ((parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size())) {
144  std::cout << "Error: parameters and errors have different number of values" << std::endl;
145  exit(1);
146  }
147 
148  fillValueError();
149 
150  sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
151  sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
152  sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
153 
154  sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
155  sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
156  sigmaPtVsEtaMinusErr_ =
157  new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
158 
159  sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
160  sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
161  sigmaMassVsPtMinusErr_ =
162  new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
163 
164  sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
165  sigmaMassVsEtaPlusErr_ =
166  new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
167  sigmaMassVsEtaMinusErr_ =
168  new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
169 
170  sigmaMassOverMassVsPt_ =
171  new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
172  sigmaMassOverMassVsPtPlusErr_ =
173  new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
174  sigmaMassOverMassVsPtMinusErr_ =
175  new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
176 
177  sigmaMassOverMassVsEta_ =
178  new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
179  sigmaMassOverMassVsEtaPlusErr_ =
180  new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
181  sigmaMassOverMassVsEtaMinusErr_ = new TProfile(
182  "sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
183 
184  sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
185  sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
186 }
187 
188 void ErrorsPropagationAnalyzer::fillValueError() {
189  valuePlusError_.resize(parameters_.size());
190  valueMinusError_.resize(parameters_.size());
191 
192  std::vector<double>::const_iterator parIt = parameters_.begin();
193  std::vector<double>::const_iterator errIt = errors_.begin();
194  std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
195  int i = 0;
196  for (; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i) {
197  valuePlusError_[i] = *parIt + (*errIt) * (*errFactorIt);
198  valueMinusError_[i] = *parIt - (*errIt) * (*errFactorIt);
199  }
200 }
201 
202 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer() {
203  gROOT->SetStyle("Plain");
204 
205  fillHistograms();
206 
207  TFile* outputFile = new TFile(outputFileName_, "RECREATE");
208 
209  drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
210  drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");
211 
212  drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
213  drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");
214 
215  drawHistograms(sigmaMassOverMassVsEta_,
216  sigmaMassOverMassVsEtaPlusErr_,
217  sigmaMassOverMassVsEtaMinusErr_,
218  "sigmaMassOverMassVsEta",
219  "#sigma(M)/M");
220  drawHistograms(sigmaMassOverMassVsPt_,
221  sigmaMassOverMassVsPtPlusErr_,
222  sigmaMassOverMassVsPtMinusErr_,
223  "sigmaMassOverMassVsPt",
224  "#sigma(M)/M");
225 
226  sigmaPtVsPtDiff_->Write();
227  sigmaPtVsEtaDiff_->Write();
228 
229  outputFile->Write();
230  outputFile->Close();
231 }
232 
233 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile* histo,
234  const TProfile* histoPlusErr,
235  const TProfile* histoMinusErr,
236  const TString& type,
237  const TString& yLabel) {
238  TH1D* sigmaPtVsEtaTH1D =
239  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
240  TH1D* sigmaPtVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
241  type + "PlusErr",
242  histo->GetNbinsX(),
243  histo->GetXaxis()->GetXmin(),
244  histo->GetXaxis()->GetXmax());
245  TH1D* sigmaPtVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
246  type + "MinusErr",
247  histo->GetNbinsX(),
248  histo->GetXaxis()->GetXmin(),
249  histo->GetXaxis()->GetXmax());
250 
251  TH1D* sigmaMassVsEtaTH1D =
252  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
253  TH1D* sigmaMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
254  type + "PlusErr",
255  histo->GetNbinsX(),
256  histo->GetXaxis()->GetXmin(),
257  histo->GetXaxis()->GetXmax());
258  TH1D* sigmaMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
259  type + "MinusErr",
260  histo->GetNbinsX(),
261  histo->GetXaxis()->GetXmin(),
262  histo->GetXaxis()->GetXmax());
263 
264  TH1D* sigmaMassOverMassVsEtaTH1D =
265  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
266  TH1D* sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
267  type + "PlusErr",
268  histo->GetNbinsX(),
269  histo->GetXaxis()->GetXmin(),
270  histo->GetXaxis()->GetXmax());
271  TH1D* sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
272  type + "MinusErr",
273  histo->GetNbinsX(),
274  histo->GetXaxis()->GetXmin(),
275  histo->GetXaxis()->GetXmax());
276 
277  TCanvas* canvas = new TCanvas("canvas" + type, "canvas" + type, 1000, 800);
278  for (int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin) {
279  sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
280  sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
281  sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
282  }
283  int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
284  // Draw TGraph with asymmetric errors
285  Double_t* values = sigmaPtVsEtaTH1D->GetArray();
286  Double_t* valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
287  Double_t* valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
288  double* posErrors = new double[numBins];
289  double* negErrors = new double[numBins];
290 
291  TGraphAsymmErrors* graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
292  TGraph* graph = new TGraph(sigmaPtVsEtaTH1D);
293 
294  for (int i = 1; i <= numBins; ++i) {
295  // std::cout << "filling " << i << std::endl;
296  posErrors[i - 1] = valuesPlus[i] - values[i];
297  if (valuesMinus[i] < 0)
298  negErrors[i - 1] = values[i];
299  else
300  negErrors[i - 1] = values[i] - valuesMinus[i];
301 
302  graphAsymmErrors->SetPointEYlow(i - 1, negErrors[i - 1]);
303  graphAsymmErrors->SetPointEYhigh(i - 1, posErrors[i - 1]);
304  }
305 
306  canvas->Draw();
307 
308  delete[] posErrors;
309  delete[] negErrors;
310 
311  graphAsymmErrors->SetTitle("");
312  graphAsymmErrors->SetFillColor(kGray);
313  graphAsymmErrors->Draw("A2");
314  TString title(type);
315  if (type == "Eta")
316  title = "#eta";
317  graphAsymmErrors->GetXaxis()->SetTitle(title);
318  graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
319  graph->Draw();
320 
321  // graph->SetLineColor(kGray);
322  // graph->SetMarkerColor(kBlack);
323  // graph->Draw("AP");
324 
325  // if( debug_ ) {
326  // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
327  // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
328  // }
329  // else {
330  // sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
331  // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
332  // sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
333  // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
334  // }
335  // sigmaPtVsEtaPlusErrTH1D->Draw();
336  // sigmaPtVsEtaTH1D->Draw("SAME");
337  // sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
338 
339  sigmaPtVsEtaPlusErrTH1D->Write();
340  sigmaPtVsEtaTH1D->Write();
341  sigmaPtVsEtaMinusErrTH1D->Write();
342 
343  sigmaPtVsEtaPlusErr_->Write();
344  sigmaPtVsEta_->Write();
345  sigmaPtVsEtaMinusErr_->Write();
346 
347  // Mass
348  sigmaMassVsEtaPlusErrTH1D->Write();
349  sigmaMassVsEtaTH1D->Write();
350  sigmaMassVsEtaMinusErrTH1D->Write();
351 
352  sigmaMassOverMassVsEtaPlusErrTH1D->Write();
353  sigmaMassOverMassVsEtaTH1D->Write();
354  sigmaMassOverMassVsEtaMinusErrTH1D->Write();
355 
356  sigmaMassVsEtaPlusErr_->Write();
357  sigmaMassVsEta_->Write();
358  sigmaMassVsEtaMinusErr_->Write();
359 
360  sigmaMassOverMassVsEtaPlusErr_->Write();
361  sigmaMassOverMassVsEta_->Write();
362  sigmaMassOverMassVsEtaMinusErr_->Write();
363 
364  canvas->Write();
365 }
366 
367 // ------------ method called to for each event ------------
368 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {}
369 
370 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
371  const lorentzVector& mu2,
372  const std::vector<double>& parval,
373  const double& sigmaPt1,
374  const double& sigmaPt2) {
375  double* p = new double[(int)(parval.size())];
376  std::vector<double>::const_iterator it = parval.begin();
377  int id = 0;
378  for (; it != parval.end(); ++it, ++id) {
379  p[id] = *it;
380  }
381  double massRes = massResolution(mu1, mu2, p, sigmaPt1, sigmaPt2);
382  delete[] p;
383  return massRes;
384 }
385 
386 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
387  const lorentzVector& mu2,
388  double* parval,
389  const double& sigmaPt1,
390  const double& sigmaPt2) {
391  double mass = (mu1 + mu2).mass();
392  double pt1 = mu1.Pt();
393  double phi1 = mu1.Phi();
394  double eta1 = mu1.Eta();
395  double theta1 = 2 * atan(exp(-eta1));
396  double pt2 = mu2.Pt();
397  double phi2 = mu2.Phi();
398  double eta2 = mu2.Eta();
399  double theta2 = 2 * atan(exp(-eta2));
400 
401  // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
402  // -----------------------------------------------------------------------------------------------------------
403  double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
404  sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
405  (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
406  pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
407  mass;
408  double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
409  sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
410  (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
411  pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
412  mass;
413  double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
414  double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
415  double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
416  sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
417  (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
418  pt1 * pt2 * cos(theta2) / sin(theta2)) /
419  mass;
420  double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
421  sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
422  (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
423  pt2 * pt1 * cos(theta1) / sin(theta1)) /
424  mass;
425 
426  // Resolution parameters:
427  // ----------------------
428  // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
429  // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
430  double sigma_pt1 = sigmaPt1;
431  double sigma_pt2 = sigmaPt2;
432  double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt1, eta1, parval);
433  double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt2, eta2, parval);
434  double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt1, eta1, parval);
435  double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt2, eta2, parval);
436  double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2(pt1, eta1, pt2, eta2, parval);
437 
438  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
439  // multiply it by pt.
440  double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
441  std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
442  std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2) +
443  2 * dmdpt1 * dmdpt2 * cov_pt1pt2);
444 
445  return mass_res;
446 }
447 
448 void ErrorsPropagationAnalyzer::fillHistograms() {
449  std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
450  RootTreeHandler rootTreeHandler;
451 
452  typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
453  MuonPairVector savedPair;
454  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
455  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
456  // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
457 
458  resolutionFunctionBase<std::vector<double> >* resolutionFunctionForVec = resolutionFunctionVecService(resolFitType_);
459  MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType_);
460  MuScleFitUtils::debugMassResol_ = false;
461  MuScleFitUtils::debug = 0;
462  MuScleFitUtils::resfind = std::vector<int>(6, 0);
463 
464  SigmaPt sigmaPt(parameters_, errors_);
465  SigmaPtDiff sigmaPtDiff;
466 
467  // Loop on all the pairs
468  unsigned int i = 0;
469  MuonPairVector::iterator it = savedPair.begin();
470  std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
471  for (; it != savedPair.end(); ++it, ++i) {
472  double pt1 = it->first.pt();
473  double eta1 = it->first.eta();
474  double pt2 = it->second.pt();
475  double eta2 = it->second.eta();
476 
477  if (debug_) {
478  std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
479  }
480  // double fabsEta1 = fabs(eta1);
481  // double fabsEta2 = fabs(eta2);
482 
483  if (pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0)
484  continue;
485 
486  // double sigmaPt1 = sigmaPt( eta1 );
487  // double sigmaPt2 = sigmaPt( eta2 );
488 
489  // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
490  // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
491  // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
492  // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
493 
494  double sigmaPt1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, parameters_);
495  double sigmaPt2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, parameters_);
496  double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
497  double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
498  double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
499  double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
500 
501  // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
502  // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
503  // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
504  double sigmaMass = massResolution(it->first, it->second, parameters_, sigmaPt1, sigmaPt2);
505  double sigmaMassPlusErr = massResolution(it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2);
506  double sigmaMassMinusErr = massResolution(it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2);
507 
508  double mass = (it->first + it->second).mass();
509 
510  if (debug_) {
511  std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
512  std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
513  std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
514  }
515 
516  // Protections from nans
517  if (pt1 != pt1)
518  continue;
519  if (pt2 != pt2)
520  continue;
521  if (eta1 != eta1)
522  continue;
523  if (eta2 != eta2)
524  continue;
525  if (sigmaPt1 != sigmaPt1)
526  continue;
527  if (sigmaPt2 != sigmaPt2)
528  continue;
529  if (sigmaPtPlusErr1 != sigmaPtPlusErr1)
530  continue;
531  if (sigmaPtPlusErr2 != sigmaPtPlusErr2)
532  continue;
533  if (sigmaPtMinusErr1 != sigmaPtMinusErr1)
534  continue;
535  if (sigmaPtMinusErr2 != sigmaPtMinusErr2)
536  continue;
537  if (sigmaMass != sigmaMass)
538  continue;
539  if (sigmaMassPlusErr != sigmaMassPlusErr)
540  continue;
541  if (sigmaMassMinusErr != sigmaMassMinusErr)
542  continue;
543  if (mass == 0)
544  continue;
545 
546  std::cout << "Muon pair number " << i << std::endl;
547 
548  // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
549  // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
550  // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
551 
552  if ((pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_)) {
553  sigmaPtVsPt_->Fill(pt1, sigmaPt1);
554  sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
555  sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
556  sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
557  sigmaMassVsPt_->Fill(pt1, sigmaMass * mass);
558  sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr * mass);
559  sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr * mass);
560  sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
561  sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
562  sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
563 
564  sigmaPtVsEta_->Fill(eta1, sigmaPt1);
565  sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
566  sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
567  sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
568  sigmaMassVsEta_->Fill(eta1, sigmaMass * mass);
569  sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr * mass);
570  sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr * mass);
571  sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
572  sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
573  sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
574  }
575  if ((pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_)) {
576  sigmaPtVsPt_->Fill(pt2, sigmaPt2);
577  sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
578  sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
579  sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
580  sigmaMassVsPt_->Fill(pt2, sigmaMass * mass);
581  sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr * mass);
582  sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr * mass);
583  sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
584  sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
585  sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
586 
587  sigmaPtVsEta_->Fill(eta2, sigmaPt2);
588  sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
589  sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
590  sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
591  sigmaMassVsEta_->Fill(eta2, sigmaMass * mass);
592  sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr * mass);
593  sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr * mass);
594  sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
595  sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
596  sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
597  }
598  }
599 }
600 
601 //define this as a plug-in
602 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
download_sqlite_cfg.outputFile
tuple outputFile
Definition: download_sqlite_cfg.py:5
mps_fire.i
i
Definition: mps_fire.py:428
SigmaPtDiff
Returns ( sigmaPt/Pt(data)^2 - sigmaPt/Pt(MC)^2 )
Definition: SigmaPtDiff.h:61
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:21
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ErrorsPropagationAnalyzer::valueMinusError_
std::vector< double > valueMinusError_
Definition: ErrorsPropagationAnalyzer.cc:92
MuScleFitUtils::debugMassResol_
static bool debugMassResol_
Definition: MuScleFitUtils.h:268
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MuScleFitUtils::debug
static int debug
Definition: MuScleFitUtils.h:130
HLT_FULL_cff.eta2
tuple eta2
Definition: HLT_FULL_cff.py:14111
RootTreeHandler::readTree
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)
Definition: RootTreeHandler.h:86
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
etaMax_
etaMax_(conf.getParameter< double >("etaMax"))
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
ErrorsPropagationAnalyzer::valuePlusError_
std::vector< double > valuePlusError_
Definition: ErrorsPropagationAnalyzer.cc:91
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
svgfig.canvas
def canvas
Definition: svgfig.py:482
resolutionFunctionService
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:70
resolutionFunctionBase::sigmaPt
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
ptMin_
ptMin_(conf.getParameter< double >("ptMin"))
iEvent
int iEvent
Definition: GenABIO.cc:224
resolutionFunctionVecService
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
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
resolutionFunctionBase::covPt1Pt2
virtual double covPt1Pt2(const double &pt1, const double &eta1, const double &pt2, const double &eta2, const T &parval)
Definition: Functions.h:834
SigmaPtDiff::squaredDiff
double squaredDiff(const double &eta)
Definition: SigmaPtDiff.h:149
ErrorsPropagationAnalyzer::massResolution
double massResolution(const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval, const double &sigmaPt1, const double &sigmaPt2)
Modified method to take into account the error.
Definition: ErrorsPropagationAnalyzer.cc:370
resolutionFunctionBase::sigmaCotgTh
virtual double sigmaCotgTh(const double &pt, const double &eta, const T &parval)=0
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
ErrorsPropagationAnalyzer::drawHistograms
void drawHistograms(const TProfile *histo, const TProfile *histoPlusErr, const TProfile *histoMinusErr, const TString &type, const TString &yLabel)
Definition: ErrorsPropagationAnalyzer.cc:233
MuonPairVector
std::vector< std::pair< lorentzVector, lorentzVector > > MuonPairVector
Definition: RootTreeHandler.h:13
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
MuScleFitUtils::mMu2
static const double mMu2
Definition: MuScleFitUtils.h:139
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer
ErrorsPropagationAnalyzer(const edm::ParameterSet &)
Definition: ErrorsPropagationAnalyzer.cc:123
resolutionFunctionBase::sigmaPhi
virtual double sigmaPhi(const double &pt, const double &eta, const T &parval)=0
MuScleFitUtils::resfind
static std::vector< int > resfind
Definition: MuScleFitUtils.h:210
HLT_FULL_cff.eta1
tuple eta1
Definition: HLT_FULL_cff.py:14110
resolutionFunctionBase::sigmaPtError
virtual double sigmaPtError(const double &pt, const double &eta, const T &parval, const T &parError)
Definition: Functions.h:831
ErrorsPropagationAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: ErrorsPropagationAnalyzer.cc:368
etaMin_
etaMin_(conf.getParameter< double >("etaMin"))
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
MuScleFitUtils::resolutionFunction
static resolutionFunctionBase< double * > * resolutionFunction
Definition: MuScleFitUtils.h:153
Generated for CMSSW Reference Manual by doxygen 1.8.5