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 #ifndef ERRORSANALYZER_CC
2 #define ERRORSANALYZER_CC
3 
4 #include "ErrorsPropagationAnalyzer.h"
5 
6 ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig)
7  : treeFileName_(iConfig.getParameter<std::string>("InputFileName")),
8  resolFitType_(iConfig.getParameter<int>("ResolFitType")),
9  maxEvents_(iConfig.getParameter<int>("MaxEvents")),
10  outputFileName_(iConfig.getParameter<std::string>("OutputFileName")),
11  ptBins_(iConfig.getParameter<int>("PtBins")),
12  ptMin_(iConfig.getParameter<double>("PtMin")),
13  ptMax_(iConfig.getParameter<double>("PtMax")),
14  etaBins_(iConfig.getParameter<int>("EtaBins")),
15  etaMin_(iConfig.getParameter<double>("EtaMin")),
16  etaMax_(iConfig.getParameter<double>("EtaMax")),
17  debug_(iConfig.getParameter<bool>("Debug")),
18  ptMinCut_(iConfig.getUntrackedParameter<double>("PtMinCut", 0.)),
19  ptMaxCut_(iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.)),
20  etaMinCut_(iConfig.getUntrackedParameter<double>("EtaMinCut", 0.)),
21  etaMaxCut_(iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.)) {
22  parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
23  errors_ = iConfig.getParameter<std::vector<double> >("Errors");
24  errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
25 
26  if ((parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size())) {
27  std::cout << "Error: parameters and errors have different number of values" << std::endl;
28  exit(1);
29  }
30 
31  fillValueError();
32 
33  sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
34  sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
35  sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
36 
37  sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
38  sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
39  sigmaPtVsEtaMinusErr_ =
40  new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
41 
42  sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
43  sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
44  sigmaMassVsPtMinusErr_ =
45  new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
46 
47  sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
48  sigmaMassVsEtaPlusErr_ =
49  new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
50  sigmaMassVsEtaMinusErr_ =
51  new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
52 
53  sigmaMassOverMassVsPt_ =
54  new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
55  sigmaMassOverMassVsPtPlusErr_ =
56  new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
57  sigmaMassOverMassVsPtMinusErr_ =
58  new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
59 
60  sigmaMassOverMassVsEta_ =
61  new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
62  sigmaMassOverMassVsEtaPlusErr_ =
63  new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
64  sigmaMassOverMassVsEtaMinusErr_ = new TProfile(
65  "sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
66 
67  sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
68  sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
69 }
70 
71 void ErrorsPropagationAnalyzer::fillValueError() {
72  valuePlusError_.resize(parameters_.size());
73  valueMinusError_.resize(parameters_.size());
74 
75  std::vector<double>::const_iterator parIt = parameters_.begin();
76  std::vector<double>::const_iterator errIt = errors_.begin();
77  std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
78  int i = 0;
79  for (; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i) {
80  valuePlusError_[i] = *parIt + (*errIt) * (*errFactorIt);
81  valueMinusError_[i] = *parIt - (*errIt) * (*errFactorIt);
82  }
83 }
84 
85 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer() {
86  gROOT->SetStyle("Plain");
87 
88  fillHistograms();
89 
90  TFile* outputFile = new TFile(outputFileName_, "RECREATE");
91 
92  drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
93  drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");
94 
95  drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
96  drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");
97 
98  drawHistograms(sigmaMassOverMassVsEta_,
99  sigmaMassOverMassVsEtaPlusErr_,
100  sigmaMassOverMassVsEtaMinusErr_,
101  "sigmaMassOverMassVsEta",
102  "#sigma(M)/M");
103  drawHistograms(sigmaMassOverMassVsPt_,
104  sigmaMassOverMassVsPtPlusErr_,
105  sigmaMassOverMassVsPtMinusErr_,
106  "sigmaMassOverMassVsPt",
107  "#sigma(M)/M");
108 
109  sigmaPtVsPtDiff_->Write();
110  sigmaPtVsEtaDiff_->Write();
111 
112  outputFile->Write();
113  outputFile->Close();
114 }
115 
116 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile* histo,
117  const TProfile* histoPlusErr,
118  const TProfile* histoMinusErr,
119  const TString& type,
120  const TString& yLabel) {
121  TH1D* sigmaPtVsEtaTH1D =
122  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
123  TH1D* sigmaPtVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
124  type + "PlusErr",
125  histo->GetNbinsX(),
126  histo->GetXaxis()->GetXmin(),
127  histo->GetXaxis()->GetXmax());
128  TH1D* sigmaPtVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
129  type + "MinusErr",
130  histo->GetNbinsX(),
131  histo->GetXaxis()->GetXmin(),
132  histo->GetXaxis()->GetXmax());
133 
134  TH1D* sigmaMassVsEtaTH1D =
135  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
136  TH1D* sigmaMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
137  type + "PlusErr",
138  histo->GetNbinsX(),
139  histo->GetXaxis()->GetXmin(),
140  histo->GetXaxis()->GetXmax());
141  TH1D* sigmaMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
142  type + "MinusErr",
143  histo->GetNbinsX(),
144  histo->GetXaxis()->GetXmin(),
145  histo->GetXaxis()->GetXmax());
146 
147  TH1D* sigmaMassOverMassVsEtaTH1D =
148  new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
149  TH1D* sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
150  type + "PlusErr",
151  histo->GetNbinsX(),
152  histo->GetXaxis()->GetXmin(),
153  histo->GetXaxis()->GetXmax());
154  TH1D* sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
155  type + "MinusErr",
156  histo->GetNbinsX(),
157  histo->GetXaxis()->GetXmin(),
158  histo->GetXaxis()->GetXmax());
159 
160  TCanvas* canvas = new TCanvas("canvas" + type, "canvas" + type, 1000, 800);
161  for (int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin) {
162  sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
163  sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
164  sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
165  }
166  int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
167  // Draw TGraph with asymmetric errors
168  Double_t* values = sigmaPtVsEtaTH1D->GetArray();
169  Double_t* valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
170  Double_t* valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
171  double* posErrors = new double[numBins];
172  double* negErrors = new double[numBins];
173 
174  TGraphAsymmErrors* graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
175  TGraph* graph = new TGraph(sigmaPtVsEtaTH1D);
176 
177  for (int i = 1; i <= numBins; ++i) {
178  // std::cout << "filling " << i << std::endl;
179  posErrors[i - 1] = valuesPlus[i] - values[i];
180  if (valuesMinus[i] < 0)
181  negErrors[i - 1] = values[i];
182  else
183  negErrors[i - 1] = values[i] - valuesMinus[i];
184 
185  graphAsymmErrors->SetPointEYlow(i - 1, negErrors[i - 1]);
186  graphAsymmErrors->SetPointEYhigh(i - 1, posErrors[i - 1]);
187  }
188 
189  canvas->Draw();
190 
191  delete[] posErrors;
192  delete[] negErrors;
193 
194  graphAsymmErrors->SetTitle("");
195  graphAsymmErrors->SetFillColor(kGray);
196  graphAsymmErrors->Draw("A2");
197  TString title(type);
198  if (type == "Eta")
199  title = "#eta";
200  graphAsymmErrors->GetXaxis()->SetTitle(title);
201  graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
202  graph->Draw();
203 
204  // graph->SetLineColor(kGray);
205  // graph->SetMarkerColor(kBlack);
206  // graph->Draw("AP");
207 
208  // if( debug_ ) {
209  // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
210  // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
211  // }
212  // else {
213  // sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
214  // sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
215  // sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
216  // sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
217  // }
218  // sigmaPtVsEtaPlusErrTH1D->Draw();
219  // sigmaPtVsEtaTH1D->Draw("SAME");
220  // sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
221 
222  sigmaPtVsEtaPlusErrTH1D->Write();
223  sigmaPtVsEtaTH1D->Write();
224  sigmaPtVsEtaMinusErrTH1D->Write();
225 
226  sigmaPtVsEtaPlusErr_->Write();
227  sigmaPtVsEta_->Write();
228  sigmaPtVsEtaMinusErr_->Write();
229 
230  // Mass
231  sigmaMassVsEtaPlusErrTH1D->Write();
232  sigmaMassVsEtaTH1D->Write();
233  sigmaMassVsEtaMinusErrTH1D->Write();
234 
235  sigmaMassOverMassVsEtaPlusErrTH1D->Write();
236  sigmaMassOverMassVsEtaTH1D->Write();
237  sigmaMassOverMassVsEtaMinusErrTH1D->Write();
238 
239  sigmaMassVsEtaPlusErr_->Write();
240  sigmaMassVsEta_->Write();
241  sigmaMassVsEtaMinusErr_->Write();
242 
243  sigmaMassOverMassVsEtaPlusErr_->Write();
244  sigmaMassOverMassVsEta_->Write();
245  sigmaMassOverMassVsEtaMinusErr_->Write();
246 
247  canvas->Write();
248 }
249 
250 // ------------ method called to for each event ------------
251 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {}
252 
253 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
254  const lorentzVector& mu2,
255  const std::vector<double>& parval,
256  const double& sigmaPt1,
257  const double& sigmaPt2) {
258  double* p = new double[(int)(parval.size())];
259  std::vector<double>::const_iterator it = parval.begin();
260  int id = 0;
261  for (; it != parval.end(); ++it, ++id) {
262  p[id] = *it;
263  }
264  double massRes = massResolution(mu1, mu2, p, sigmaPt1, sigmaPt2);
265  delete[] p;
266  return massRes;
267 }
268 
269 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
270  const lorentzVector& mu2,
271  double* parval,
272  const double& sigmaPt1,
273  const double& sigmaPt2) {
274  double mass = (mu1 + mu2).mass();
275  double pt1 = mu1.Pt();
276  double phi1 = mu1.Phi();
277  double eta1 = mu1.Eta();
278  double theta1 = 2 * atan(exp(-eta1));
279  double pt2 = mu2.Pt();
280  double phi2 = mu2.Phi();
281  double eta2 = mu2.Eta();
282  double theta2 = 2 * atan(exp(-eta2));
283 
284  // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
285  // -----------------------------------------------------------------------------------------------------------
286  double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
287  sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
288  (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
289  pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
290  mass;
291  double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
292  sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
293  (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
294  pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
295  mass;
296  double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
297  double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
298  double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
299  sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
300  (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
301  pt1 * pt2 * cos(theta2) / sin(theta2)) /
302  mass;
303  double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
304  sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
305  (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
306  pt2 * pt1 * cos(theta1) / sin(theta1)) /
307  mass;
308 
309  // Resolution parameters:
310  // ----------------------
311  // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
312  // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
313  double sigma_pt1 = sigmaPt1;
314  double sigma_pt2 = sigmaPt2;
315  double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt1, eta1, parval);
316  double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt2, eta2, parval);
317  double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt1, eta1, parval);
318  double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt2, eta2, parval);
319  double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2(pt1, eta1, pt2, eta2, parval);
320 
321  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
322  // multiply it by pt.
323  double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
324  std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
325  std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2) +
326  2 * dmdpt1 * dmdpt2 * cov_pt1pt2);
327 
328  return mass_res;
329 }
330 
331 void ErrorsPropagationAnalyzer::fillHistograms() {
332  std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
333  RootTreeHandler rootTreeHandler;
334 
335  typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
336  MuonPairVector savedPair;
337  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
338  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
339  // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
340 
341  resolutionFunctionBase<std::vector<double> >* resolutionFunctionForVec = resolutionFunctionVecService(resolFitType_);
342  MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType_);
343  MuScleFitUtils::debugMassResol_ = false;
344  MuScleFitUtils::debug = 0;
345  MuScleFitUtils::resfind = std::vector<int>(6, 0);
346 
347  SigmaPt sigmaPt(parameters_, errors_);
348  SigmaPtDiff sigmaPtDiff;
349 
350  // Loop on all the pairs
351  unsigned int i = 0;
352  MuonPairVector::iterator it = savedPair.begin();
353  std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
354  for (; it != savedPair.end(); ++it, ++i) {
355  double pt1 = it->first.pt();
356  double eta1 = it->first.eta();
357  double pt2 = it->second.pt();
358  double eta2 = it->second.eta();
359 
360  if (debug_) {
361  std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
362  }
363  // double fabsEta1 = fabs(eta1);
364  // double fabsEta2 = fabs(eta2);
365 
366  if (pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0)
367  continue;
368 
369  // double sigmaPt1 = sigmaPt( eta1 );
370  // double sigmaPt2 = sigmaPt( eta2 );
371 
372  // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
373  // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
374  // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
375  // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
376 
377  double sigmaPt1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, parameters_);
378  double sigmaPt2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, parameters_);
379  double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
380  double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
381  double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
382  double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
383 
384  // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
385  // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
386  // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
387  double sigmaMass = massResolution(it->first, it->second, parameters_, sigmaPt1, sigmaPt2);
388  double sigmaMassPlusErr = massResolution(it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2);
389  double sigmaMassMinusErr = massResolution(it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2);
390 
391  double mass = (it->first + it->second).mass();
392 
393  if (debug_) {
394  std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
395  std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
396  std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
397  }
398 
399  // Protections from nans
400  if (pt1 != pt1)
401  continue;
402  if (pt2 != pt2)
403  continue;
404  if (eta1 != eta1)
405  continue;
406  if (eta2 != eta2)
407  continue;
408  if (sigmaPt1 != sigmaPt1)
409  continue;
410  if (sigmaPt2 != sigmaPt2)
411  continue;
412  if (sigmaPtPlusErr1 != sigmaPtPlusErr1)
413  continue;
414  if (sigmaPtPlusErr2 != sigmaPtPlusErr2)
415  continue;
416  if (sigmaPtMinusErr1 != sigmaPtMinusErr1)
417  continue;
418  if (sigmaPtMinusErr2 != sigmaPtMinusErr2)
419  continue;
420  if (sigmaMass != sigmaMass)
421  continue;
422  if (sigmaMassPlusErr != sigmaMassPlusErr)
423  continue;
424  if (sigmaMassMinusErr != sigmaMassMinusErr)
425  continue;
426  if (mass == 0)
427  continue;
428 
429  std::cout << "Muon pair number " << i << std::endl;
430 
431  // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
432  // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
433  // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
434 
435  if ((pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_)) {
436  sigmaPtVsPt_->Fill(pt1, sigmaPt1);
437  sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
438  sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
439  sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
440  sigmaMassVsPt_->Fill(pt1, sigmaMass * mass);
441  sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr * mass);
442  sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr * mass);
443  sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
444  sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
445  sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
446 
447  sigmaPtVsEta_->Fill(eta1, sigmaPt1);
448  sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
449  sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
450  sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
451  sigmaMassVsEta_->Fill(eta1, sigmaMass * mass);
452  sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr * mass);
453  sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr * mass);
454  sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
455  sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
456  sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
457  }
458  if ((pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_)) {
459  sigmaPtVsPt_->Fill(pt2, sigmaPt2);
460  sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
461  sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
462  sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
463  sigmaMassVsPt_->Fill(pt2, sigmaMass * mass);
464  sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr * mass);
465  sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr * mass);
466  sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
467  sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
468  sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
469 
470  sigmaPtVsEta_->Fill(eta2, sigmaPt2);
471  sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
472  sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
473  sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
474  sigmaMassVsEta_->Fill(eta2, sigmaMass * mass);
475  sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr * mass);
476  sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr * mass);
477  sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
478  sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
479  sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
480  }
481  }
482 }
483 
484 //define this as a plug-in
485 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
486 
487 #endif
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:20
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ErrorsPropagationAnalyzer::valueMinusError_
std::vector< double > valueMinusError_
Definition: ErrorsPropagationAnalyzer.h:95
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:9780
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.h:94
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:253
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:116
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:6
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:9779
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:251
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