CMS 3D CMS Logo

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