CMS 3D CMS Logo
EcalZmassClient.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: ZfitterAnalyzer
4 // Class: ZfitterAnalyzer
5 //
14 //
15 // Original Author: Vieri Candelise
16 // Created: Mon Jun 13 09:49:08 CEST 2011
17 //
18 //
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
24 #include "DQM/Physics/src/EwkDQM.h"
25 #include "DQMServices/Core/interface/DQMEDHarvester.h"
26 #include "DQMServices/Core/interface/MonitorElement.h"
27 #include "FWCore/Framework/interface/MakerMacros.h"
28 #include "FWCore/ParameterSet/interface/ParameterSet.h"
29 #include "FWCore/Utilities/interface/Exception.h"
30 #include "TH1.h"
31 #include "TMath.h"
32 #include <cmath>
33 #include <iostream>
34 #include <string>
35 
36 Double_t mybw(Double_t *, Double_t *);
37 Double_t mygauss(Double_t *, Double_t *);
38 
39 class EcalZmassClient : public DQMEDHarvester {
40 public:
41  explicit EcalZmassClient(const edm::ParameterSet &);
42  ~EcalZmassClient() override;
43 
44 private:
45  void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override;
46 
47  // ----------member data ---------------------------
48 
49  std::string prefixME_;
50 };
51 
52 EcalZmassClient::EcalZmassClient(const edm::ParameterSet &iConfig)
53  : prefixME_(iConfig.getUntrackedParameter<std::string>("prefixME", "")) {}
54 
55 EcalZmassClient::~EcalZmassClient() {}
56 
57 void EcalZmassClient::dqmEndJob(DQMStore::IBooker &_ibooker, DQMStore::IGetter &_igetter) {
58  MonitorElement *h_fitres1;
59  MonitorElement *h_fitres1bis;
60  MonitorElement *h_fitres1Chi2;
61  MonitorElement *h_fitres2;
62  MonitorElement *h_fitres2bis;
63  MonitorElement *h_fitres2Chi2;
64  MonitorElement *h_fitres3;
65  MonitorElement *h_fitres3bis;
66  MonitorElement *h_fitres3Chi2;
67 
68  _ibooker.setCurrentFolder(prefixME_ + "/Zmass");
69  h_fitres1 = _ibooker.book1D("Gaussian mean WP80 EB-EB", "Gaussian mean WP80 EB-EB", 1, 0, 1);
70  h_fitres1bis = _ibooker.book1D("Gaussian sigma WP80 EB-EB", "Gaussian sigma WP80 EB-EB", 1, 0, 1);
71  h_fitres1Chi2 = _ibooker.book1D(
72  "Gaussian Chi2 result over NDF WP80 EB-EB", "Gaussian Chi2 result over NDF WP80 EB-EB", 1, 0, 1);
73 
74  h_fitres3 = _ibooker.book1D("Gaussian mean WP80 EB-EE", "Gaussian mean result WP80 EB-EE", 1, 0, 1);
75  h_fitres3bis = _ibooker.book1D("Gaussian sigma WP80 EB-EE", "Gaussian sigma WP80 EB-EE", 1, 0, 1);
76  h_fitres3Chi2 =
77  _ibooker.book1D("Gaussian Chi2 result over NDF WP80 EB-EE", "Gaussian Chi2 result over NDF WP80 EB-EE", 1, 0, 1);
78 
79  h_fitres2 = _ibooker.book1D("Gaussian mean WP80 EE-EE", "Gaussian mean WP80 EE-EE", 1, 0, 1);
80  h_fitres2bis = _ibooker.book1D("Gaussian sigma WP80 EE-EE", "Gaussian sigma WP80 EE-EE", 1, 0, 1);
81  h_fitres2Chi2 =
82  _ibooker.book1D("Gaussian Chi2 result over NDF WP80 EE-EE", "Gaussian Chi2 result over NDF WP80 EE-EE", 1, 0, 1);
83 
84  LogTrace("EwkAnalyzer") << "Parameters initialization";
85 
86  MonitorElement *me1 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EB-EB");
87  MonitorElement *me2 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EE-EE");
88  MonitorElement *me3 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EB-EE");
89 
90  if (me1 != nullptr) {
91  TH1F *B = me1->getTH1F();
92  TH1F *R1 = h_fitres1->getTH1F();
93  int division = B->GetNbinsX();
94  float massMIN = B->GetBinLowEdge(1);
95  float massMAX = B->GetBinLowEdge(division + 1);
96  // float BIN_SIZE = B->GetBinWidth(1);
97 
98  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
99  func->SetParameter(0, 1.0);
100  func->SetParName(0, "const");
101  func->SetParameter(1, 95.0);
102  func->SetParName(1, "mean");
103  func->SetParameter(2, 5.0);
104  func->SetParName(2, "sigma");
105 
106  double stats[4];
107  R1->GetStats(stats);
108  float N = 0;
109  float mean = 0;
110  float sigma = 0;
111  N = B->GetEntries();
112 
113  try {
114  if (N != 0) {
115  B->Fit("mygauss", "QR");
116  mean = std::abs(func->GetParameter(1));
117  sigma = std::abs(func->GetParError(1));
118  }
119 
120  if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
121  N = 1;
122  mean = 0;
123  sigma = 0;
124  }
125 
126  } catch (cms::Exception &e) {
127  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
128  N = 1;
129  mean = 40;
130  sigma = 0;
131  }
132 
133  stats[0] = N;
134  stats[1] = N;
135  stats[2] = mean * N;
136  stats[3] = sigma * sigma * N + mean * mean * N;
137 
138  R1->SetEntries(N);
139  R1->PutStats(stats);
140  }
141  /*******************************************************/
142  if (me1 != nullptr) {
143  TH1F *Bbis = me1->getTH1F();
144  TH1F *R1bis = h_fitres1bis->getTH1F();
145  int division = Bbis->GetNbinsX();
146  float massMIN = Bbis->GetBinLowEdge(1);
147  float massMAX = Bbis->GetBinLowEdge(division + 1);
148  // float BIN_SIZE = B->GetBinWidth(1);
149 
150  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
151  func->SetParameter(0, 1.0);
152  func->SetParName(0, "const");
153  func->SetParameter(1, 95.0);
154  func->SetParName(1, "mean");
155  func->SetParameter(2, 5.0);
156  func->SetParName(2, "sigma");
157 
158  double stats[4];
159  R1bis->GetStats(stats);
160  float N = 0;
161  float rms = 0;
162  float rmsErr = 0;
163  N = Bbis->GetEntries();
164 
165  try {
166  if (N != 0) {
167  Bbis->Fit("mygauss", "QR");
168  rms = std::abs(func->GetParameter(2));
169  rmsErr = std::abs(func->GetParError(2));
170  }
171 
172  if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
173  N = 1;
174  rms = 0;
175  rmsErr = 0;
176  }
177 
178  } catch (cms::Exception &e) {
179  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
180  N = 1;
181  rms = 40;
182  rmsErr = 0;
183  }
184 
185  stats[0] = N;
186  stats[1] = N;
187  stats[2] = rms * N;
188  stats[3] = rmsErr * rmsErr * N + rms * rms * N;
189 
190  R1bis->SetEntries(N);
191  R1bis->PutStats(stats);
192  }
193  /****************************************/
194 
195  if (me2 != nullptr) {
196  TH1F *E = me2->getTH1F();
197  TH1F *R2 = h_fitres2->getTH1F();
198  int division = E->GetNbinsX();
199  float massMIN = E->GetBinLowEdge(1);
200  float massMAX = E->GetBinLowEdge(division + 1);
201  // float BIN_SIZE = E->GetBinWidth(1);
202 
203  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
204  func->SetParameter(0, 1.0);
205  func->SetParName(0, "const");
206  func->SetParameter(1, 95.0);
207  func->SetParName(1, "mean");
208  func->SetParameter(2, 5.0);
209  func->SetParName(2, "sigma");
210 
211  double stats[4];
212  R2->GetStats(stats);
213  float N = 0;
214  float mean = 0;
215  float sigma = 0;
216  N = E->GetEntries();
217 
218  try {
219  if (N != 0) {
220  E->Fit("mygauss", "QR");
221  mean = std::abs(func->GetParameter(1));
222  sigma = std::abs(func->GetParError(1));
223  }
224 
225  if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
226  N = 1;
227  mean = 0;
228  sigma = 0;
229  }
230 
231  } catch (cms::Exception &e) {
232  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
233  N = 1;
234  mean = 40;
235  sigma = 0;
236  }
237 
238  stats[0] = N;
239  stats[1] = N;
240  stats[2] = mean * N;
241  stats[3] = sigma * sigma * N + mean * mean * N;
242 
243  R2->SetEntries(N);
244  R2->PutStats(stats);
245  }
246  /**************************************************************************/
247 
248  if (me2 != nullptr) {
249  TH1F *Ebis = me2->getTH1F();
250  TH1F *R2bis = h_fitres2bis->getTH1F();
251  int division = Ebis->GetNbinsX();
252  float massMIN = Ebis->GetBinLowEdge(1);
253  float massMAX = Ebis->GetBinLowEdge(division + 1);
254  // float BIN_SIZE = B->GetBinWidth(1);
255 
256  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
257  func->SetParameter(0, 1.0);
258  func->SetParName(0, "const");
259  func->SetParameter(1, 95.0);
260  func->SetParName(1, "mean");
261  func->SetParameter(2, 5.0);
262  func->SetParName(2, "sigma");
263 
264  double stats[4];
265  R2bis->GetStats(stats);
266  float N = 0;
267  float rms = 0;
268  float rmsErr = 0;
269  N = Ebis->GetEntries();
270 
271  try {
272  if (N != 0) {
273  Ebis->Fit("mygauss", "QR");
274  rms = std::abs(func->GetParameter(2));
275  rmsErr = std::abs(func->GetParError(2));
276  }
277 
278  if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
279  N = 1;
280  rms = 0;
281  rmsErr = 0;
282  }
283 
284  } catch (cms::Exception &e) {
285  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
286  N = 1;
287  rms = 40;
288  rmsErr = 0;
289  }
290 
291  stats[0] = N;
292  stats[1] = N;
293  stats[2] = rms * N;
294  stats[3] = rmsErr * rmsErr * N + rms * rms * N;
295 
296  R2bis->SetEntries(N);
297  R2bis->PutStats(stats);
298  }
299  /*********************************************************************************************/
300 
301  if (me3 != nullptr) {
302  TH1F *R3 = h_fitres3->getTH1F();
303  TH1F *M = me3->getTH1F();
304  int division = M->GetNbinsX();
305  float massMIN = M->GetBinLowEdge(1);
306  float massMAX = M->GetBinLowEdge(division + 1);
307  // float BIN_SIZE = M->GetBinWidth(1);
308 
309  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
310  func->SetParameter(0, 1.0);
311  func->SetParName(0, "const");
312  func->SetParameter(1, 95.0);
313  func->SetParName(1, "mean");
314  func->SetParameter(2, 5.0);
315  func->SetParName(2, "sigma");
316 
317  double stats[4];
318  R3->GetStats(stats);
319  float N = 0;
320  float mean = 0;
321  float sigma = 0;
322  N = M->GetEntries();
323 
324  try {
325  if (N != 0) {
326  M->Fit("mygauss", "QR");
327  mean = std::abs(func->GetParameter(1));
328  sigma = std::abs(func->GetParError(1));
329  }
330  if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
331  N = 1;
332  mean = 0;
333  sigma = 0;
334  }
335 
336  } catch (cms::Exception &e) {
337  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
338  N = 1;
339  mean = 40;
340  sigma = 0;
341  }
342 
343  stats[0] = N;
344  stats[1] = N;
345  stats[2] = mean * N;
346  stats[3] = sigma * sigma * N + mean * mean * N;
347 
348  R3->SetEntries(N);
349  R3->PutStats(stats);
350  }
351  /********************************************************************************/
352 
353  if (me3 != nullptr) {
354  TH1F *Mbis = me3->getTH1F();
355  TH1F *R3bis = h_fitres3bis->getTH1F();
356  int division = Mbis->GetNbinsX();
357  float massMIN = Mbis->GetBinLowEdge(1);
358  float massMAX = Mbis->GetBinLowEdge(division + 1);
359  // float BIN_SIZE = B->GetBinWidth(1);
360 
361  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
362  func->SetParameter(0, 1.0);
363  func->SetParName(0, "const");
364  func->SetParameter(1, 95.0);
365  func->SetParName(1, "mean");
366  func->SetParameter(2, 5.0);
367  func->SetParName(2, "sigma");
368 
369  double stats[4];
370  R3bis->GetStats(stats);
371  float N = 0;
372  float rms = 0;
373  float rmsErr = 0;
374  N = Mbis->GetEntries();
375 
376  try {
377  if (N != 0) {
378  Mbis->Fit("mygauss", "QR");
379  rms = std::abs(func->GetParameter(2));
380  rmsErr = std::abs(func->GetParError(2));
381  }
382 
383  if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
384  N = 1;
385  rms = 0;
386  rmsErr = 0;
387  }
388 
389  } catch (cms::Exception &e) {
390  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
391  N = 1;
392  rms = 40;
393  rmsErr = 0;
394  }
395 
396  stats[0] = N;
397  stats[1] = N;
398  stats[2] = rms * N;
399  stats[3] = rmsErr * rmsErr * N + rms * rms * N;
400 
401  R3bis->SetEntries(N);
402  R3bis->PutStats(stats);
403  }
404 
405  /*Chi2 */
406 
407  if (me1 != nullptr) {
408  TH1F *C1 = me1->getTH1F();
409  TH1F *S1 = h_fitres1Chi2->getTH1F();
410  int division = C1->GetNbinsX();
411  float massMIN = C1->GetBinLowEdge(1);
412  float massMAX = C1->GetBinLowEdge(division + 1);
413  // float BIN_SIZE = B->GetBinWidth(1);
414 
415  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
416  func->SetParameter(0, 1.0);
417  func->SetParName(0, "const");
418  func->SetParameter(1, 95.0);
419  func->SetParName(1, "mean");
420  func->SetParameter(2, 5.0);
421  func->SetParName(2, "sigma");
422 
423  double stats[4];
424  S1->GetStats(stats);
425  float N = 0;
426  float Chi2 = 0;
427  float NDF = 0;
428  N = C1->GetEntries();
429 
430  try {
431  if (N != 0) {
432  C1->Fit("mygauss", "QR");
433  if ((func->GetNDF() != 0)) {
434  Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
435  NDF = 0.1;
436  }
437  }
438 
439  if (N == 0 || Chi2 < 0 || NDF < 0) {
440  N = 1;
441  Chi2 = 0;
442  NDF = 0;
443  }
444 
445  } catch (cms::Exception &e) {
446  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
447  N = 1;
448  Chi2 = 40;
449  NDF = 0;
450  }
451 
452  stats[0] = N;
453  stats[1] = N;
454  stats[2] = Chi2 * N;
455  stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
456 
457  S1->SetEntries(N);
458  S1->PutStats(stats);
459  }
460  /**********************************************/
461 
462  if (me2 != nullptr) {
463  TH1F *C2 = me2->getTH1F();
464  TH1F *S2 = h_fitres2Chi2->getTH1F();
465  int division = C2->GetNbinsX();
466  float massMIN = C2->GetBinLowEdge(1);
467  float massMAX = C2->GetBinLowEdge(division + 1);
468  // float BIN_SIZE = B->GetBinWidth(1);
469 
470  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
471  func->SetParameter(0, 1.0);
472  func->SetParName(0, "const");
473  func->SetParameter(1, 95.0);
474  func->SetParName(1, "mean");
475  func->SetParameter(2, 5.0);
476  func->SetParName(2, "sigma");
477 
478  double stats[4];
479  S2->GetStats(stats);
480  float N = 0;
481  float Chi2 = 0;
482  float NDF = 0;
483  N = C2->GetEntries();
484 
485  try {
486  if (N != 0) {
487  C2->Fit("mygauss", "QR");
488  if (func->GetNDF() != 0) {
489  Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
490  NDF = 0.1;
491  }
492  }
493 
494  if (N == 0 || Chi2 < 0 || NDF < 0) {
495  N = 1;
496  Chi2 = 0;
497  NDF = 0;
498  }
499 
500  } catch (cms::Exception &e) {
501  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
502  N = 1;
503  Chi2 = 40;
504  NDF = 0;
505  }
506 
507  stats[0] = N;
508  stats[1] = N;
509  stats[2] = Chi2 * N;
510  stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
511 
512  S2->SetEntries(N);
513  S2->PutStats(stats);
514  }
515  /**************************************************************************/
516  if (me3 != nullptr) {
517  TH1F *C3 = me3->getTH1F();
518  TH1F *S3 = h_fitres3Chi2->getTH1F();
519  int division = C3->GetNbinsX();
520  float massMIN = C3->GetBinLowEdge(1);
521  float massMAX = C3->GetBinLowEdge(division + 1);
522  // float BIN_SIZE = B->GetBinWidth(1);
523 
524  TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
525  func->SetParameter(0, 1.0);
526  func->SetParName(0, "const");
527  func->SetParameter(1, 95.0);
528  func->SetParName(1, "mean");
529  func->SetParameter(2, 5.0);
530  func->SetParName(2, "sigma");
531 
532  double stats[4];
533  S3->GetStats(stats);
534  float N = 0;
535  float Chi2 = 0;
536  float NDF = 0;
537  N = C3->GetEntries();
538 
539  try {
540  if (N != 0) {
541  C3->Fit("mygauss", "QR");
542  if ((func->GetNDF() != 0)) {
543  Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
544  NDF = 0.1;
545  }
546  }
547 
548  if (N == 0 || Chi2 < 0 || NDF < 0) {
549  N = 1;
550  Chi2 = 0;
551  NDF = 0;
552  }
553 
554  } catch (cms::Exception &e) {
555  edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
556  N = 1;
557  Chi2 = 40;
558  NDF = 0;
559  }
560 
561  stats[0] = N;
562  stats[1] = N;
563  stats[2] = Chi2 * N;
564  stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
565 
566  S3->SetEntries(N);
567  S3->PutStats(stats);
568  }
569 }
570 
571 // Breit-Wigner function
572 Double_t mybw(Double_t *x, Double_t *par) {
573  Double_t arg1 = 14.0 / 22.0; // 2 over pi
574  Double_t arg2 = par[1] * par[1] * par[2] * par[2]; // Gamma=par[2] M=par[1]
575  Double_t arg3 = ((x[0] * x[0]) - (par[1] * par[1])) * ((x[0] * x[0]) - (par[1] * par[1]));
576  Double_t arg4 = x[0] * x[0] * x[0] * x[0] * ((par[2] * par[2]) / (par[1] * par[1]));
577  return par[0] * arg1 * arg2 / (arg3 + arg4);
578 }
579 
580 // Gaussian
581 Double_t mygauss(Double_t *x, Double_t *par) {
582  Double_t arg = 0;
583  if (par[2] < 0)
584  par[2] = -par[2]; // par[2]: sigma
585  if (par[2] != 0)
586  arg = (x[0] - par[1]) / par[2]; // par[1]: mean
587  return par[0] * TMath::Exp(-0.5 * arg * arg) / (TMath::Sqrt(2 * TMath::Pi()) * par[2]); // par[0] is constant
588 }
589 
590 DEFINE_FWK_MODULE(EcalZmassClient);
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
MonitorElement::getTH1F
TH1F * getTH1F() const
Definition: MonitorElement.cc:1483
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
EcalZmassClient::~EcalZmassClient
~EcalZmassClient() override
Definition: EcalZmassClient.cc:55
jets_cff.func
func
Definition: jets_cff.py:495
cms::Exception::what
char const * what() const override
Definition: Exception.cc:103
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
EcalZmassClient::dqmEndJob
void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override
Definition: EcalZmassClient.cc:57
funct::arg
A arg
Definition: Factorize.h:38
mybw
Double_t mybw(Double_t *, Double_t *)
Definition: EcalZmassClient.cc:572
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:268
DQMStore::IBooker::book1D
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:106
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:671
DQMStore::IGetter::get
MonitorElement * get(std::string const &path)
Definition: DQMStore.cc:303
N
#define N
Definition: blowfish.cc:9
mygauss
Double_t mygauss(Double_t *, Double_t *)
Definition: EcalZmassClient.cc:581
Chi2
Definition: Chi2.h:15
EcalZmassClient::EcalZmassClient
EcalZmassClient(const edm::ParameterSet &)
Definition: EcalZmassClient.cc:52
EcalZmassClient::prefixME_
std::string prefixME_
Definition: EcalZmassClient.cc:49
Generated for CMSSW Reference Manual by doxygen 1.8.11