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