Ecal/plugins/EcalZmassClient.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ZfitterAnalyzer
0004 // Class:      ZfitterAnalyzer
0005 //
0006 /**\class ZfitterAnalyzer ZfitterAnalyzer.cc
0007  Zfitter/ZfitterAnalyzer/src/ZfitterAnalyzer.cc
0008
0009  Description: [one line class summary]
0010
0011  Implementation:
0012      [Notes on implementation]
0013 */
0014 //
0015 // Original Author:  Vieri Candelise
0016 //         Created:  Mon Jun 13 09:49:08 CEST 2011
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022
0023 // user include files
0024 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0025 #include "DQMServices/Core/interface/DQMStore.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/Utilities/interface/Exception.h"
0029 #include "TH1.h"
0030 #include "TMath.h"
0031 #include <cmath>
0032 #include <iostream>
0033 #include <string>
0034
0035 Double_t mybw(Double_t *, Double_t *);
0036 Double_t mygauss(Double_t *, Double_t *);
0037
0038 class EcalZmassClient : public DQMEDHarvester {
0039 public:
0040   explicit EcalZmassClient(const edm::ParameterSet &);
0041   ~EcalZmassClient() override;
0042
0043 private:
0044   void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override;
0045
0046   // ----------member data ---------------------------
0047
0048   std::string prefixME_;
0049 };
0050
0051 EcalZmassClient::EcalZmassClient(const edm::ParameterSet &iConfig)
0052     : prefixME_(iConfig.getUntrackedParameter<std::string>("prefixME", "")) {}
0053
0054 EcalZmassClient::~EcalZmassClient() {}
0055
0056 void EcalZmassClient::dqmEndJob(DQMStore::IBooker &_ibooker, DQMStore::IGetter &_igetter) {
0057   MonitorElement *h_fitres1;
0058   MonitorElement *h_fitres1bis;
0059   MonitorElement *h_fitres1Chi2;
0060   MonitorElement *h_fitres2;
0061   MonitorElement *h_fitres2bis;
0062   MonitorElement *h_fitres2Chi2;
0063   MonitorElement *h_fitres3;
0064   MonitorElement *h_fitres3bis;
0065   MonitorElement *h_fitres3Chi2;
0066
0067   _ibooker.setCurrentFolder(prefixME_ + "/Zmass");
0068   h_fitres1 = _ibooker.book1D("Gaussian mean WP80 EB-EB", "Gaussian mean WP80 EB-EB", 1, 0, 1);
0069   h_fitres1bis = _ibooker.book1D("Gaussian sigma WP80 EB-EB", "Gaussian sigma WP80 EB-EB", 1, 0, 1);
0070   h_fitres1Chi2 = _ibooker.book1D(
0071       "Gaussian Chi2 result over NDF  WP80 EB-EB", "Gaussian Chi2 result over NDF  WP80 EB-EB", 1, 0, 1);
0072
0073   h_fitres3 = _ibooker.book1D("Gaussian mean WP80 EB-EE", "Gaussian mean result WP80 EB-EE", 1, 0, 1);
0074   h_fitres3bis = _ibooker.book1D("Gaussian sigma WP80 EB-EE", "Gaussian sigma WP80 EB-EE", 1, 0, 1);
0075   h_fitres3Chi2 =
0076       _ibooker.book1D("Gaussian Chi2 result over NDF WP80 EB-EE", "Gaussian Chi2 result over NDF WP80 EB-EE", 1, 0, 1);
0077
0078   h_fitres2 = _ibooker.book1D("Gaussian mean WP80 EE-EE", "Gaussian mean WP80 EE-EE", 1, 0, 1);
0079   h_fitres2bis = _ibooker.book1D("Gaussian sigma WP80 EE-EE", "Gaussian sigma WP80 EE-EE", 1, 0, 1);
0080   h_fitres2Chi2 =
0081       _ibooker.book1D("Gaussian Chi2 result over NDF WP80 EE-EE", "Gaussian Chi2 result over NDF WP80 EE-EE", 1, 0, 1);
0082
0083   LogTrace("EwkAnalyzer") << "Parameters initialization";
0084
0085   MonitorElement *me1 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EB-EB");
0086   MonitorElement *me2 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EE-EE");
0087   MonitorElement *me3 = _igetter.get(prefixME_ + "/Zmass/Z peak - WP80 EB-EE");
0088
0089   if (me1 != nullptr) {
0090     TH1F *B = me1->getTH1F();
0091     TH1F *R1 = h_fitres1->getTH1F();
0092     int division = B->GetNbinsX();
0093     float massMIN = B->GetBinLowEdge(1);
0094     float massMAX = B->GetBinLowEdge(division + 1);
0095     // float BIN_SIZE = B->GetBinWidth(1);
0096
0097     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0098     func->SetParameter(0, 1.0);
0099     func->SetParName(0, "const");
0100     func->SetParameter(1, 95.0);
0101     func->SetParName(1, "mean");
0102     func->SetParameter(2, 5.0);
0103     func->SetParName(2, "sigma");
0104
0105     double stats[4];
0106     R1->GetStats(stats);
0107     float N = 0;
0108     float mean = 0;
0109     float sigma = 0;
0110     N = B->GetEntries();
0111
0112     try {
0113       if (N != 0) {
0114         B->Fit("mygauss", "QR");
0115         mean = std::abs(func->GetParameter(1));
0116         sigma = std::abs(func->GetParError(1));
0117       }
0118
0119       if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
0120         N = 1;
0121         mean = 0;
0122         sigma = 0;
0123       }
0124
0125     } catch (cms::Exception &e) {
0126       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0127       N = 1;
0128       mean = 40;
0129       sigma = 0;
0130     }
0131
0132     stats[0] = N;
0133     stats[1] = N;
0134     stats[2] = mean * N;
0135     stats[3] = sigma * sigma * N + mean * mean * N;
0136
0137     R1->SetEntries(N);
0138     R1->PutStats(stats);
0139   }
0140   /*******************************************************/
0141   if (me1 != nullptr) {
0142     TH1F *Bbis = me1->getTH1F();
0143     TH1F *R1bis = h_fitres1bis->getTH1F();
0144     int division = Bbis->GetNbinsX();
0145     float massMIN = Bbis->GetBinLowEdge(1);
0146     float massMAX = Bbis->GetBinLowEdge(division + 1);
0147     // float BIN_SIZE = B->GetBinWidth(1);
0148
0149     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0150     func->SetParameter(0, 1.0);
0151     func->SetParName(0, "const");
0152     func->SetParameter(1, 95.0);
0153     func->SetParName(1, "mean");
0154     func->SetParameter(2, 5.0);
0155     func->SetParName(2, "sigma");
0156
0157     double stats[4];
0158     R1bis->GetStats(stats);
0159     float N = 0;
0160     float rms = 0;
0161     float rmsErr = 0;
0162     N = Bbis->GetEntries();
0163
0164     try {
0165       if (N != 0) {
0166         Bbis->Fit("mygauss", "QR");
0167         rms = std::abs(func->GetParameter(2));
0168         rmsErr = std::abs(func->GetParError(2));
0169       }
0170
0171       if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
0172         N = 1;
0173         rms = 0;
0174         rmsErr = 0;
0175       }
0176
0177     } catch (cms::Exception &e) {
0178       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0179       N = 1;
0180       rms = 40;
0181       rmsErr = 0;
0182     }
0183
0184     stats[0] = N;
0185     stats[1] = N;
0186     stats[2] = rms * N;
0187     stats[3] = rmsErr * rmsErr * N + rms * rms * N;
0188
0189     R1bis->SetEntries(N);
0190     R1bis->PutStats(stats);
0191   }
0192   /****************************************/
0193
0194   if (me2 != nullptr) {
0195     TH1F *E = me2->getTH1F();
0196     TH1F *R2 = h_fitres2->getTH1F();
0197     int division = E->GetNbinsX();
0198     float massMIN = E->GetBinLowEdge(1);
0199     float massMAX = E->GetBinLowEdge(division + 1);
0200     // float BIN_SIZE = E->GetBinWidth(1);
0201
0202     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0203     func->SetParameter(0, 1.0);
0204     func->SetParName(0, "const");
0205     func->SetParameter(1, 95.0);
0206     func->SetParName(1, "mean");
0207     func->SetParameter(2, 5.0);
0208     func->SetParName(2, "sigma");
0209
0210     double stats[4];
0211     R2->GetStats(stats);
0212     float N = 0;
0213     float mean = 0;
0214     float sigma = 0;
0215     N = E->GetEntries();
0216
0217     try {
0218       if (N != 0) {
0219         E->Fit("mygauss", "QR");
0220         mean = std::abs(func->GetParameter(1));
0221         sigma = std::abs(func->GetParError(1));
0222       }
0223
0224       if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
0225         N = 1;
0226         mean = 0;
0227         sigma = 0;
0228       }
0229
0230     } catch (cms::Exception &e) {
0231       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0232       N = 1;
0233       mean = 40;
0234       sigma = 0;
0235     }
0236
0237     stats[0] = N;
0238     stats[1] = N;
0239     stats[2] = mean * N;
0240     stats[3] = sigma * sigma * N + mean * mean * N;
0241
0242     R2->SetEntries(N);
0243     R2->PutStats(stats);
0244   }
0245   /**************************************************************************/
0246
0247   if (me2 != nullptr) {
0248     TH1F *Ebis = me2->getTH1F();
0249     TH1F *R2bis = h_fitres2bis->getTH1F();
0250     int division = Ebis->GetNbinsX();
0251     float massMIN = Ebis->GetBinLowEdge(1);
0252     float massMAX = Ebis->GetBinLowEdge(division + 1);
0253     // float BIN_SIZE = B->GetBinWidth(1);
0254
0255     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0256     func->SetParameter(0, 1.0);
0257     func->SetParName(0, "const");
0258     func->SetParameter(1, 95.0);
0259     func->SetParName(1, "mean");
0260     func->SetParameter(2, 5.0);
0261     func->SetParName(2, "sigma");
0262
0263     double stats[4];
0264     R2bis->GetStats(stats);
0265     float N = 0;
0266     float rms = 0;
0267     float rmsErr = 0;
0268     N = Ebis->GetEntries();
0269
0270     try {
0271       if (N != 0) {
0272         Ebis->Fit("mygauss", "QR");
0273         rms = std::abs(func->GetParameter(2));
0274         rmsErr = std::abs(func->GetParError(2));
0275       }
0276
0277       if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
0278         N = 1;
0279         rms = 0;
0280         rmsErr = 0;
0281       }
0282
0283     } catch (cms::Exception &e) {
0284       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0285       N = 1;
0286       rms = 40;
0287       rmsErr = 0;
0288     }
0289
0290     stats[0] = N;
0291     stats[1] = N;
0292     stats[2] = rms * N;
0293     stats[3] = rmsErr * rmsErr * N + rms * rms * N;
0294
0295     R2bis->SetEntries(N);
0296     R2bis->PutStats(stats);
0297   }
0298   /*********************************************************************************************/
0299
0300   if (me3 != nullptr) {
0301     TH1F *R3 = h_fitres3->getTH1F();
0302     TH1F *M = me3->getTH1F();
0303     int division = M->GetNbinsX();
0304     float massMIN = M->GetBinLowEdge(1);
0305     float massMAX = M->GetBinLowEdge(division + 1);
0306     // float BIN_SIZE = M->GetBinWidth(1);
0307
0308     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0309     func->SetParameter(0, 1.0);
0310     func->SetParName(0, "const");
0311     func->SetParameter(1, 95.0);
0312     func->SetParName(1, "mean");
0313     func->SetParameter(2, 5.0);
0314     func->SetParName(2, "sigma");
0315
0316     double stats[4];
0317     R3->GetStats(stats);
0318     float N = 0;
0319     float mean = 0;
0320     float sigma = 0;
0321     N = M->GetEntries();
0322
0323     try {
0324       if (N != 0) {
0325         M->Fit("mygauss", "QR");
0326         mean = std::abs(func->GetParameter(1));
0327         sigma = std::abs(func->GetParError(1));
0328       }
0329       if (N == 0 || mean < 50 || mean > 100 || sigma <= 0 || sigma > 20) {
0330         N = 1;
0331         mean = 0;
0332         sigma = 0;
0333       }
0334
0335     } catch (cms::Exception &e) {
0336       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0337       N = 1;
0338       mean = 40;
0339       sigma = 0;
0340     }
0341
0342     stats[0] = N;
0343     stats[1] = N;
0344     stats[2] = mean * N;
0345     stats[3] = sigma * sigma * N + mean * mean * N;
0346
0347     R3->SetEntries(N);
0348     R3->PutStats(stats);
0349   }
0350   /********************************************************************************/
0351
0352   if (me3 != nullptr) {
0353     TH1F *Mbis = me3->getTH1F();
0354     TH1F *R3bis = h_fitres3bis->getTH1F();
0355     int division = Mbis->GetNbinsX();
0356     float massMIN = Mbis->GetBinLowEdge(1);
0357     float massMAX = Mbis->GetBinLowEdge(division + 1);
0358     // float BIN_SIZE = B->GetBinWidth(1);
0359
0360     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0361     func->SetParameter(0, 1.0);
0362     func->SetParName(0, "const");
0363     func->SetParameter(1, 95.0);
0364     func->SetParName(1, "mean");
0365     func->SetParameter(2, 5.0);
0366     func->SetParName(2, "sigma");
0367
0368     double stats[4];
0369     R3bis->GetStats(stats);
0370     float N = 0;
0371     float rms = 0;
0372     float rmsErr = 0;
0373     N = Mbis->GetEntries();
0374
0375     try {
0376       if (N != 0) {
0377         Mbis->Fit("mygauss", "QR");
0378         rms = std::abs(func->GetParameter(2));
0379         rmsErr = std::abs(func->GetParError(2));
0380       }
0381
0382       if (N == 0 || rms < 0 || rms > 50 || rmsErr <= 0 || rmsErr > 50) {
0383         N = 1;
0384         rms = 0;
0385         rmsErr = 0;
0386       }
0387
0388     } catch (cms::Exception &e) {
0389       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0390       N = 1;
0391       rms = 40;
0392       rmsErr = 0;
0393     }
0394
0395     stats[0] = N;
0396     stats[1] = N;
0397     stats[2] = rms * N;
0398     stats[3] = rmsErr * rmsErr * N + rms * rms * N;
0399
0400     R3bis->SetEntries(N);
0401     R3bis->PutStats(stats);
0402   }
0403
0404   /*Chi2 */
0405
0406   if (me1 != nullptr) {
0407     TH1F *C1 = me1->getTH1F();
0408     TH1F *S1 = h_fitres1Chi2->getTH1F();
0409     int division = C1->GetNbinsX();
0410     float massMIN = C1->GetBinLowEdge(1);
0411     float massMAX = C1->GetBinLowEdge(division + 1);
0412     // float BIN_SIZE = B->GetBinWidth(1);
0413
0414     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0415     func->SetParameter(0, 1.0);
0416     func->SetParName(0, "const");
0417     func->SetParameter(1, 95.0);
0418     func->SetParName(1, "mean");
0419     func->SetParameter(2, 5.0);
0420     func->SetParName(2, "sigma");
0421
0422     double stats[4];
0423     S1->GetStats(stats);
0424     float N = 0;
0425     float Chi2 = 0;
0426     float NDF = 0;
0427     N = C1->GetEntries();
0428
0429     try {
0430       if (N != 0) {
0431         C1->Fit("mygauss", "QR");
0432         if ((func->GetNDF() != 0)) {
0433           Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
0434           NDF = 0.1;
0435         }
0436       }
0437
0438       if (N == 0 || Chi2 < 0 || NDF < 0) {
0439         N = 1;
0440         Chi2 = 0;
0441         NDF = 0;
0442       }
0443
0444     } catch (cms::Exception &e) {
0445       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0446       N = 1;
0447       Chi2 = 40;
0448       NDF = 0;
0449     }
0450
0451     stats[0] = N;
0452     stats[1] = N;
0453     stats[2] = Chi2 * N;
0454     stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
0455
0456     S1->SetEntries(N);
0457     S1->PutStats(stats);
0458   }
0459   /**********************************************/
0460
0461   if (me2 != nullptr) {
0462     TH1F *C2 = me2->getTH1F();
0463     TH1F *S2 = h_fitres2Chi2->getTH1F();
0464     int division = C2->GetNbinsX();
0465     float massMIN = C2->GetBinLowEdge(1);
0466     float massMAX = C2->GetBinLowEdge(division + 1);
0467     // float BIN_SIZE = B->GetBinWidth(1);
0468
0469     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0470     func->SetParameter(0, 1.0);
0471     func->SetParName(0, "const");
0472     func->SetParameter(1, 95.0);
0473     func->SetParName(1, "mean");
0474     func->SetParameter(2, 5.0);
0475     func->SetParName(2, "sigma");
0476
0477     double stats[4];
0478     S2->GetStats(stats);
0479     float N = 0;
0480     float Chi2 = 0;
0481     float NDF = 0;
0482     N = C2->GetEntries();
0483
0484     try {
0485       if (N != 0) {
0486         C2->Fit("mygauss", "QR");
0487         if (func->GetNDF() != 0) {
0488           Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
0489           NDF = 0.1;
0490         }
0491       }
0492
0493       if (N == 0 || Chi2 < 0 || NDF < 0) {
0494         N = 1;
0495         Chi2 = 0;
0496         NDF = 0;
0497       }
0498
0499     } catch (cms::Exception &e) {
0500       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0501       N = 1;
0502       Chi2 = 40;
0503       NDF = 0;
0504     }
0505
0506     stats[0] = N;
0507     stats[1] = N;
0508     stats[2] = Chi2 * N;
0509     stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
0510
0511     S2->SetEntries(N);
0512     S2->PutStats(stats);
0513   }
0514   /**************************************************************************/
0515   if (me3 != nullptr) {
0516     TH1F *C3 = me3->getTH1F();
0517     TH1F *S3 = h_fitres3Chi2->getTH1F();
0518     int division = C3->GetNbinsX();
0519     float massMIN = C3->GetBinLowEdge(1);
0520     float massMAX = C3->GetBinLowEdge(division + 1);
0521     // float BIN_SIZE = B->GetBinWidth(1);
0522
0523     TF1 *func = new TF1("mygauss", mygauss, massMIN, massMAX, 3);
0524     func->SetParameter(0, 1.0);
0525     func->SetParName(0, "const");
0526     func->SetParameter(1, 95.0);
0527     func->SetParName(1, "mean");
0528     func->SetParameter(2, 5.0);
0529     func->SetParName(2, "sigma");
0530
0531     double stats[4];
0532     S3->GetStats(stats);
0533     float N = 0;
0534     float Chi2 = 0;
0535     float NDF = 0;
0536     N = C3->GetEntries();
0537
0538     try {
0539       if (N != 0) {
0540         C3->Fit("mygauss", "QR");
0541         if ((func->GetNDF() != 0)) {
0542           Chi2 = std::abs(func->GetChisquare()) / std::abs(func->GetNDF());
0543           NDF = 0.1;
0544         }
0545       }
0546
0547       if (N == 0 || Chi2 < 0 || NDF < 0) {
0548         N = 1;
0549         Chi2 = 0;
0550         NDF = 0;
0551       }
0552
0553     } catch (cms::Exception &e) {
0554       edm::LogError("ZFitter") << "[Zfitter]: Exception when fitting..." << e.what();
0555       N = 1;
0556       Chi2 = 40;
0557       NDF = 0;
0558     }
0559
0560     stats[0] = N;
0561     stats[1] = N;
0562     stats[2] = Chi2 * N;
0563     stats[3] = NDF * NDF * N + Chi2 * Chi2 * N;
0564
0565     S3->SetEntries(N);
0566     S3->PutStats(stats);
0567   }
0568 }
0569
0570 // Breit-Wigner function
0571 Double_t mybw(Double_t *x, Double_t *par) {
0572   Double_t arg1 = 14.0 / 22.0;                        // 2 over pi
0573   Double_t arg2 = par[1] * par[1] * par[2] * par[2];  // Gamma=par[2]  M=par[1]
0574   Double_t arg3 = ((x[0] * x[0]) - (par[1] * par[1])) * ((x[0] * x[0]) - (par[1] * par[1]));
0575   Double_t arg4 = x[0] * x[0] * x[0] * x[0] * ((par[2] * par[2]) / (par[1] * par[1]));
0576   return par[0] * arg1 * arg2 / (arg3 + arg4);
0577 }
0578
0579 // Gaussian
0580 Double_t mygauss(Double_t *x, Double_t *par) {
0581   Double_t arg = 0;
0582   if (par[2] < 0)
0583     par[2] = -par[2];  // par[2]: sigma
0584   if (par[2] != 0)
0585     arg = (x[0] - par[1]) / par[2];                                                        // par[1]: mean
0586   return par[0] * TMath::Exp(-0.5 * arg * arg) / (TMath::Sqrt(2 * TMath::Pi()) * par[2]);  // par[0] is constant
0587 }
0588
0589 DEFINE_FWK_MODULE(EcalZmassClient);