MomentumScaleCalibration/plugins/ErrorsAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ErrorsAnalyzer
0004 // Class:      ErrorsAnalyzer
0005 //
0006 /**\class ErrorsAnalyzer ErrorsAnalyzer.cc MuonAnalysis/MomentumScaleCalibration/plugins/ErrorsAnalyzer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Marco De Mattia
0015 //         Created:  Thu Sep 11 12:16:00 CEST 2008
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <vector>
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include <TH1D.h>
0031 #include <TProfile.h>
0032 #include <TString.h>
0033 #include <TCanvas.h>
0034 #include <TGraphAsymmErrors.h>
0035 #include <TROOT.h>
0036
0037 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
0038 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
0039 #include "MuScleFitUtils.h"
0040
0041 //
0042 // class declaration
0043 //
0044
0045 class ErrorsAnalyzer : public edm::one::EDAnalyzer<> {
0046 public:
0047   explicit ErrorsAnalyzer(const edm::ParameterSet&);
0048   ~ErrorsAnalyzer() override;
0049
0050 private:
0051   void analyze(const edm::Event&, const edm::EventSetup&) override;
0052   void fillHistograms();
0053   void drawHistograms(const TProfile* histo,
0054                       const TProfile* histoPlusErr,
0055                       const TProfile* histoMinusErr,
0056                       const TString& type);
0057   void fillValueError();
0058   void endJob() override {}
0059
0060   TString treeFileName_;
0061   int resolFitType_;
0062   uint32_t maxEvents_;
0063   TString outputFileName_;
0064   int ptBins_;
0065   double ptMin_;
0066   double ptMax_;
0067   int etaBins_;
0068   double etaMin_;
0069   double etaMax_;
0070   bool debug_;
0071
0072   std::vector<double> parameters_;
0073   std::vector<double> errors_;
0074   std::vector<int> errorFactors_;
0075
0076   std::vector<double> valuePlusError_;
0077   std::vector<double> valueMinusError_;
0078
0079   TProfile* sigmaPtVsEta_;
0080   TProfile* sigmaPtVsEtaPlusErr_;
0081   TProfile* sigmaPtVsEtaMinusErr_;
0082
0083   TProfile* sigmaPtVsPt_;
0084   TProfile* sigmaPtVsPtPlusErr_;
0085   TProfile* sigmaPtVsPtMinusErr_;
0086
0087   // Mass resolution
0088   TProfile* sigmaMassVsEta_;
0089   TProfile* sigmaMassVsEtaPlusErr_;
0090   TProfile* sigmaMassVsEtaMinusErr_;
0091
0092   TProfile* sigmaMassVsPt_;
0093   TProfile* sigmaMassVsPtPlusErr_;
0094   TProfile* sigmaMassVsPtMinusErr_;
0095 };
0096
0097 ErrorsAnalyzer::ErrorsAnalyzer(const edm::ParameterSet& iConfig)
0098     : treeFileName_(iConfig.getParameter<std::string>("InputFileName")),
0099       resolFitType_(iConfig.getParameter<int>("ResolFitType")),
0100       maxEvents_(iConfig.getParameter<int>("MaxEvents")),
0101       outputFileName_(iConfig.getParameter<std::string>("OutputFileName")),
0102       ptBins_(iConfig.getParameter<int>("PtBins")),
0103       ptMin_(iConfig.getParameter<double>("PtMin")),
0104       ptMax_(iConfig.getParameter<double>("PtMax")),
0105       etaBins_(iConfig.getParameter<int>("EtaBins")),
0106       etaMin_(iConfig.getParameter<double>("EtaMin")),
0107       etaMax_(iConfig.getParameter<double>("EtaMax")),
0108       debug_(iConfig.getParameter<bool>("Debug")) {
0109   parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
0110   errors_ = iConfig.getParameter<std::vector<double> >("Errors");
0111   errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
0112
0113   if ((parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size())) {
0114     std::cout << "Error: parameters and errors have different number of values" << std::endl;
0115     exit(1);
0116   }
0117
0118   fillValueError();
0119
0120   sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
0121   sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
0122   sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
0123
0124   sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
0125   sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
0126   sigmaPtVsEtaMinusErr_ =
0127       new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
0128
0129   sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
0130   sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
0131   sigmaMassVsPtMinusErr_ =
0132       new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
0133
0134   sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
0135   sigmaMassVsEtaPlusErr_ =
0136       new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
0137   sigmaMassVsEtaMinusErr_ =
0138       new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
0139 }
0140
0141 void ErrorsAnalyzer::fillValueError() {
0142   valuePlusError_.resize(parameters_.size());
0143   valueMinusError_.resize(parameters_.size());
0144
0145   std::vector<double>::const_iterator parIt = parameters_.begin();
0146   std::vector<double>::const_iterator errIt = errors_.begin();
0147   std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
0148   int i = 0;
0149   for (; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i) {
0150     valuePlusError_[i] = *parIt + (*errIt) * (*errFactorIt);
0151     valueMinusError_[i] = *parIt - (*errIt) * (*errFactorIt);
0152   }
0153 }
0154
0155 ErrorsAnalyzer::~ErrorsAnalyzer() {
0156   gROOT->SetStyle("Plain");
0157
0158   fillHistograms();
0159
0160   TFile* outputFile = new TFile(outputFileName_, "RECREATE");
0161
0162   drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta");
0163   drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt");
0164
0165   drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta");
0166   drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt");
0167
0168   outputFile->Write();
0169   outputFile->Close();
0170 }
0171
0172 void ErrorsAnalyzer::drawHistograms(const TProfile* histo,
0173                                     const TProfile* histoPlusErr,
0174                                     const TProfile* histoMinusErr,
0175                                     const TString& type) {
0176   TH1D* sigmaPtVsEtaTH1D =
0177       new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0178   TH1D* sigmaPtVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
0179                                            type + "PlusErr",
0180                                            histo->GetNbinsX(),
0181                                            histo->GetXaxis()->GetXmin(),
0182                                            histo->GetXaxis()->GetXmax());
0183   TH1D* sigmaPtVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
0184                                             type + "MinusErr",
0185                                             histo->GetNbinsX(),
0186                                             histo->GetXaxis()->GetXmin(),
0187                                             histo->GetXaxis()->GetXmax());
0188
0189   TH1D* sigmaMassVsEtaTH1D =
0190       new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0191   TH1D* sigmaMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
0192                                              type + "PlusErr",
0193                                              histo->GetNbinsX(),
0194                                              histo->GetXaxis()->GetXmin(),
0195                                              histo->GetXaxis()->GetXmax());
0196   TH1D* sigmaMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
0197                                               type + "MinusErr",
0198                                               histo->GetNbinsX(),
0199                                               histo->GetXaxis()->GetXmin(),
0200                                               histo->GetXaxis()->GetXmax());
0201
0202   TCanvas* canvas = new TCanvas("canvas" + type, "canvas" + type, 1000, 800);
0203   for (int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin) {
0204     sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
0205     sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
0206     sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
0207   }
0208   int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
0209   // Draw TGraph with asymmetric errors
0210   Double_t* values = sigmaPtVsEtaTH1D->GetArray();
0211   Double_t* valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
0212   Double_t* valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
0213   double* posErrors = new double[numBins];
0214   double* negErrors = new double[numBins];
0215
0216   TGraphAsymmErrors* graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
0217   TGraph* graph = new TGraph(sigmaPtVsEtaTH1D);
0218
0219   for (int i = 1; i <= numBins; ++i) {
0220     // std::cout << "filling " << i << std::endl;
0221     posErrors[i - 1] = valuesPlus[i] - values[i];
0222     if (valuesMinus[i] < 0)
0223       negErrors[i - 1] = values[i];
0224     else
0225       negErrors[i - 1] = values[i] - valuesMinus[i];
0226
0227     graphAsymmErrors->SetPointEYlow(i - 1, negErrors[i - 1]);
0228     graphAsymmErrors->SetPointEYhigh(i - 1, posErrors[i - 1]);
0229   }
0230
0231   canvas->Draw();
0232
0233   delete[] posErrors;
0234   delete[] negErrors;
0235
0236   graphAsymmErrors->SetTitle("");
0237   graphAsymmErrors->SetFillColor(kGray);
0238   graphAsymmErrors->Draw("A2");
0239   TString title(type);
0240   if (type == "Eta")
0241     title = "#eta";
0242   graphAsymmErrors->GetXaxis()->SetTitle(title);
0243   graphAsymmErrors->GetYaxis()->SetTitle("#sigmaPt/Pt");
0244   graph->Draw();
0245
0246   //  graph->SetLineColor(kGray);
0247   //  graph->SetMarkerColor(kBlack);
0248   //  graph->Draw("AP");
0249
0250   //   if( debug_ ) {
0251   //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
0252   //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
0253   //   }
0254   //   else {
0255   //     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
0256   //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
0257   //     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
0258   //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
0259   //   }
0260   //   sigmaPtVsEtaPlusErrTH1D->Draw();
0261   //   sigmaPtVsEtaTH1D->Draw("SAME");
0262   //   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
0263
0264   sigmaPtVsEtaPlusErrTH1D->Write();
0265   sigmaPtVsEtaTH1D->Write();
0266   sigmaPtVsEtaMinusErrTH1D->Write();
0267
0268   sigmaPtVsEtaPlusErr_->Write();
0269   sigmaPtVsEta_->Write();
0270   sigmaPtVsEtaMinusErr_->Write();
0271
0272   // Mass
0273   sigmaMassVsEtaPlusErrTH1D->Write();
0274   sigmaMassVsEtaTH1D->Write();
0275   sigmaMassVsEtaMinusErrTH1D->Write();
0276
0277   sigmaMassVsEtaPlusErr_->Write();
0278   sigmaMassVsEta_->Write();
0279   sigmaMassVsEtaMinusErr_->Write();
0280
0281   canvas->Write();
0282 }
0283
0284 // ------------ method called to for each event  ------------
0285 void ErrorsAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {}
0286
0287 void ErrorsAnalyzer::fillHistograms() {
0288   std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
0289   RootTreeHandler rootTreeHandler;
0290
0291   typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
0292   MuonPairVector savedPair;
0293   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0294   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
0295   // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
0296
0297   resolutionFunctionBase<std::vector<double> >* resolutionFunctionForVec = resolutionFunctionVecService(resolFitType_);
0298   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType_);
0299   MuScleFitUtils::debugMassResol_ = false;
0300   MuScleFitUtils::debug = 0;
0301   MuScleFitUtils::resfind = std::vector<int>(6, 0);
0302
0303   // Loop on all the pairs
0304   unsigned int i = 0;
0305   MuonPairVector::iterator it = savedPair.begin();
0306   std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
0307   for (; it != savedPair.end(); ++it, ++i) {
0308     double pt1 = it->first.pt();
0309     double eta1 = it->first.eta();
0310     double pt2 = it->second.pt();
0311     double eta2 = it->second.eta();
0312
0313     if (debug_) {
0314       std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
0315     }
0316
0317     if (pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0)
0318       continue;
0319
0320     double sigmaPt1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, parameters_);
0321     double sigmaPt2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, parameters_);
0322     double sigmaPtPlusErr1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, valuePlusError_);
0323     double sigmaPtPlusErr2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, valuePlusError_);
0324     double sigmaPtMinusErr1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, valueMinusError_);
0325     double sigmaPtMinusErr2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, valueMinusError_);
0326
0327     double sigmaMass = MuScleFitUtils::massResolution(it->first, it->second, parameters_);
0328     double sigmaMassPlusErr = MuScleFitUtils::massResolution(it->first, it->second, valuePlusError_);
0329     double sigmaMassMinusErr = MuScleFitUtils::massResolution(it->first, it->second, valueMinusError_);
0330
0331     if (debug_) {
0332       std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
0333       std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
0334       std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
0335     }
0336
0337     // Protections from nans
0338     if (pt1 != pt1)
0339       continue;
0340     if (pt2 != pt2)
0341       continue;
0342     if (eta1 != eta1)
0343       continue;
0344     if (eta2 != eta2)
0345       continue;
0346     if (sigmaPt1 != sigmaPt1)
0347       continue;
0348     if (sigmaPt2 != sigmaPt2)
0349       continue;
0350     if (sigmaPtPlusErr1 != sigmaPtPlusErr1)
0351       continue;
0352     if (sigmaPtPlusErr2 != sigmaPtPlusErr2)
0353       continue;
0354     if (sigmaPtMinusErr1 != sigmaPtMinusErr1)
0355       continue;
0356     if (sigmaPtMinusErr2 != sigmaPtMinusErr2)
0357       continue;
0358     if (sigmaMass != sigmaMass)
0359       continue;
0360     if (sigmaMassPlusErr != sigmaMassPlusErr)
0361       continue;
0362     if (sigmaMassMinusErr != sigmaMassMinusErr)
0363       continue;
0364
0365     std::cout << "Muon pair number " << i << std::endl;
0366
0367     // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
0368     // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
0369     // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
0370
0371     sigmaPtVsPt_->Fill(pt1, sigmaPt1);
0372     sigmaPtVsPt_->Fill(pt2, sigmaPt2);
0373     sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
0374     sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
0375     sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
0376     sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
0377
0378     sigmaPtVsEta_->Fill(eta1, sigmaPt1);
0379     sigmaPtVsEta_->Fill(eta2, sigmaPt2);
0380     sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
0381     sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
0382     sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
0383     sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
0384
0385     sigmaMassVsPt_->Fill(pt1, sigmaMass);
0386     sigmaMassVsPt_->Fill(pt2, sigmaMass);
0387     sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
0388     sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
0389     sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
0390     sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
0391
0392     sigmaMassVsEta_->Fill(eta1, sigmaMass);
0393     sigmaMassVsEta_->Fill(eta2, sigmaMass);
0394     sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
0395     sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
0396     sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
0397     sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
0398   }
0399 }
0400
0401 //define this as a plug-in
0402 DEFINE_FWK_MODULE(ErrorsAnalyzer);