Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:12 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ErrorsPropagationAnalyzer
0004 // Class:      ErrorsPropagationAnalyzer
0005 //
0006 /**\class ErrorsPropagationAnalyzer ErrorsPropagationAnalyzer.cc MuonAnalysis/MomentumScaleCalibration/plugins/ErrorsPropagationAnalyzer.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 #include "MuonAnalysis/MomentumScaleCalibration/interface/SigmaPtDiff.h"
0041 
0042 //
0043 // class declaration
0044 //
0045 
0046 class ErrorsPropagationAnalyzer : public edm::one::EDAnalyzer<> {
0047 public:
0048   explicit ErrorsPropagationAnalyzer(const edm::ParameterSet&);
0049   ~ErrorsPropagationAnalyzer() override;
0050 
0051 private:
0052   void analyze(const edm::Event&, const edm::EventSetup&) override;
0053   void fillHistograms();
0054   void drawHistograms(const TProfile* histo,
0055                       const TProfile* histoPlusErr,
0056                       const TProfile* histoMinusErr,
0057                       const TString& type,
0058                       const TString& yLabel);
0059   void fillValueError();
0060   void endJob() override {}
0061   /// Modified method to take into account the error
0062   double massResolution(const lorentzVector& mu1,
0063                         const lorentzVector& mu2,
0064                         const std::vector<double>& parval,
0065                         const double& sigmaPt1,
0066                         const double& sigmaPt2);
0067   double massResolution(const lorentzVector& mu1,
0068                         const lorentzVector& mu2,
0069                         double* parval,
0070                         const double& sigmaPt1,
0071                         const double& sigmaPt2);
0072 
0073   TString treeFileName_;
0074   int resolFitType_;
0075   uint32_t maxEvents_;
0076   TString outputFileName_;
0077   int ptBins_;
0078   double ptMin_;
0079   double ptMax_;
0080   int etaBins_;
0081   double etaMin_;
0082   double etaMax_;
0083   bool debug_;
0084 
0085   double ptMinCut_, ptMaxCut_, etaMinCut_, etaMaxCut_;
0086 
0087   std::vector<double> parameters_;
0088   std::vector<double> errors_;
0089   std::vector<int> errorFactors_;
0090 
0091   std::vector<double> valuePlusError_;
0092   std::vector<double> valueMinusError_;
0093 
0094   TProfile* sigmaPtVsEta_;
0095   TProfile* sigmaPtVsEtaPlusErr_;
0096   TProfile* sigmaPtVsEtaMinusErr_;
0097 
0098   TProfile* sigmaPtVsPt_;
0099   TProfile* sigmaPtVsPtPlusErr_;
0100   TProfile* sigmaPtVsPtMinusErr_;
0101 
0102   TProfile* sigmaPtVsEtaDiff_;
0103   TProfile* sigmaPtVsPtDiff_;
0104 
0105   // Mass resolution
0106   TProfile* sigmaMassVsEta_;
0107   TProfile* sigmaMassVsEtaPlusErr_;
0108   TProfile* sigmaMassVsEtaMinusErr_;
0109 
0110   TProfile* sigmaMassVsPt_;
0111   TProfile* sigmaMassVsPtPlusErr_;
0112   TProfile* sigmaMassVsPtMinusErr_;
0113 
0114   TProfile* sigmaMassOverMassVsEta_;
0115   TProfile* sigmaMassOverMassVsEtaPlusErr_;
0116   TProfile* sigmaMassOverMassVsEtaMinusErr_;
0117 
0118   TProfile* sigmaMassOverMassVsPt_;
0119   TProfile* sigmaMassOverMassVsPtPlusErr_;
0120   TProfile* sigmaMassOverMassVsPtMinusErr_;
0121 };
0122 
0123 ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig)
0124     : treeFileName_(iConfig.getParameter<std::string>("InputFileName")),
0125       resolFitType_(iConfig.getParameter<int>("ResolFitType")),
0126       maxEvents_(iConfig.getParameter<int>("MaxEvents")),
0127       outputFileName_(iConfig.getParameter<std::string>("OutputFileName")),
0128       ptBins_(iConfig.getParameter<int>("PtBins")),
0129       ptMin_(iConfig.getParameter<double>("PtMin")),
0130       ptMax_(iConfig.getParameter<double>("PtMax")),
0131       etaBins_(iConfig.getParameter<int>("EtaBins")),
0132       etaMin_(iConfig.getParameter<double>("EtaMin")),
0133       etaMax_(iConfig.getParameter<double>("EtaMax")),
0134       debug_(iConfig.getParameter<bool>("Debug")),
0135       ptMinCut_(iConfig.getUntrackedParameter<double>("PtMinCut", 0.)),
0136       ptMaxCut_(iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.)),
0137       etaMinCut_(iConfig.getUntrackedParameter<double>("EtaMinCut", 0.)),
0138       etaMaxCut_(iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.)) {
0139   parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
0140   errors_ = iConfig.getParameter<std::vector<double> >("Errors");
0141   errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
0142 
0143   if ((parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size())) {
0144     std::cout << "Error: parameters and errors have different number of values" << std::endl;
0145     exit(1);
0146   }
0147 
0148   fillValueError();
0149 
0150   sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
0151   sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
0152   sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
0153 
0154   sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
0155   sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
0156   sigmaPtVsEtaMinusErr_ =
0157       new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
0158 
0159   sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
0160   sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
0161   sigmaMassVsPtMinusErr_ =
0162       new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
0163 
0164   sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
0165   sigmaMassVsEtaPlusErr_ =
0166       new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
0167   sigmaMassVsEtaMinusErr_ =
0168       new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
0169 
0170   sigmaMassOverMassVsPt_ =
0171       new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
0172   sigmaMassOverMassVsPtPlusErr_ =
0173       new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
0174   sigmaMassOverMassVsPtMinusErr_ =
0175       new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
0176 
0177   sigmaMassOverMassVsEta_ =
0178       new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
0179   sigmaMassOverMassVsEtaPlusErr_ =
0180       new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
0181   sigmaMassOverMassVsEtaMinusErr_ = new TProfile(
0182       "sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
0183 
0184   sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
0185   sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
0186 }
0187 
0188 void ErrorsPropagationAnalyzer::fillValueError() {
0189   valuePlusError_.resize(parameters_.size());
0190   valueMinusError_.resize(parameters_.size());
0191 
0192   std::vector<double>::const_iterator parIt = parameters_.begin();
0193   std::vector<double>::const_iterator errIt = errors_.begin();
0194   std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
0195   int i = 0;
0196   for (; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i) {
0197     valuePlusError_[i] = *parIt + (*errIt) * (*errFactorIt);
0198     valueMinusError_[i] = *parIt - (*errIt) * (*errFactorIt);
0199   }
0200 }
0201 
0202 ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer() {
0203   gROOT->SetStyle("Plain");
0204 
0205   fillHistograms();
0206 
0207   TFile* outputFile = new TFile(outputFileName_, "RECREATE");
0208 
0209   drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
0210   drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");
0211 
0212   drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
0213   drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");
0214 
0215   drawHistograms(sigmaMassOverMassVsEta_,
0216                  sigmaMassOverMassVsEtaPlusErr_,
0217                  sigmaMassOverMassVsEtaMinusErr_,
0218                  "sigmaMassOverMassVsEta",
0219                  "#sigma(M)/M");
0220   drawHistograms(sigmaMassOverMassVsPt_,
0221                  sigmaMassOverMassVsPtPlusErr_,
0222                  sigmaMassOverMassVsPtMinusErr_,
0223                  "sigmaMassOverMassVsPt",
0224                  "#sigma(M)/M");
0225 
0226   sigmaPtVsPtDiff_->Write();
0227   sigmaPtVsEtaDiff_->Write();
0228 
0229   outputFile->Write();
0230   outputFile->Close();
0231 }
0232 
0233 void ErrorsPropagationAnalyzer::drawHistograms(const TProfile* histo,
0234                                                const TProfile* histoPlusErr,
0235                                                const TProfile* histoMinusErr,
0236                                                const TString& type,
0237                                                const TString& yLabel) {
0238   TH1D* sigmaPtVsEtaTH1D =
0239       new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0240   TH1D* sigmaPtVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
0241                                            type + "PlusErr",
0242                                            histo->GetNbinsX(),
0243                                            histo->GetXaxis()->GetXmin(),
0244                                            histo->GetXaxis()->GetXmax());
0245   TH1D* sigmaPtVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
0246                                             type + "MinusErr",
0247                                             histo->GetNbinsX(),
0248                                             histo->GetXaxis()->GetXmin(),
0249                                             histo->GetXaxis()->GetXmax());
0250 
0251   TH1D* sigmaMassVsEtaTH1D =
0252       new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0253   TH1D* sigmaMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
0254                                              type + "PlusErr",
0255                                              histo->GetNbinsX(),
0256                                              histo->GetXaxis()->GetXmin(),
0257                                              histo->GetXaxis()->GetXmax());
0258   TH1D* sigmaMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
0259                                               type + "MinusErr",
0260                                               histo->GetNbinsX(),
0261                                               histo->GetXaxis()->GetXmin(),
0262                                               histo->GetXaxis()->GetXmax());
0263 
0264   TH1D* sigmaMassOverMassVsEtaTH1D =
0265       new TH1D(type, type, histo->GetNbinsX(), histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0266   TH1D* sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type + "PlusErr",
0267                                                      type + "PlusErr",
0268                                                      histo->GetNbinsX(),
0269                                                      histo->GetXaxis()->GetXmin(),
0270                                                      histo->GetXaxis()->GetXmax());
0271   TH1D* sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type + "MinusErr",
0272                                                       type + "MinusErr",
0273                                                       histo->GetNbinsX(),
0274                                                       histo->GetXaxis()->GetXmin(),
0275                                                       histo->GetXaxis()->GetXmax());
0276 
0277   TCanvas* canvas = new TCanvas("canvas" + type, "canvas" + type, 1000, 800);
0278   for (int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin) {
0279     sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
0280     sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
0281     sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
0282   }
0283   int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
0284   // Draw TGraph with asymmetric errors
0285   Double_t* values = sigmaPtVsEtaTH1D->GetArray();
0286   Double_t* valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
0287   Double_t* valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
0288   double* posErrors = new double[numBins];
0289   double* negErrors = new double[numBins];
0290 
0291   TGraphAsymmErrors* graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
0292   TGraph* graph = new TGraph(sigmaPtVsEtaTH1D);
0293 
0294   for (int i = 1; i <= numBins; ++i) {
0295     // std::cout << "filling " << i << std::endl;
0296     posErrors[i - 1] = valuesPlus[i] - values[i];
0297     if (valuesMinus[i] < 0)
0298       negErrors[i - 1] = values[i];
0299     else
0300       negErrors[i - 1] = values[i] - valuesMinus[i];
0301 
0302     graphAsymmErrors->SetPointEYlow(i - 1, negErrors[i - 1]);
0303     graphAsymmErrors->SetPointEYhigh(i - 1, posErrors[i - 1]);
0304   }
0305 
0306   canvas->Draw();
0307 
0308   delete[] posErrors;
0309   delete[] negErrors;
0310 
0311   graphAsymmErrors->SetTitle("");
0312   graphAsymmErrors->SetFillColor(kGray);
0313   graphAsymmErrors->Draw("A2");
0314   TString title(type);
0315   if (type == "Eta")
0316     title = "#eta";
0317   graphAsymmErrors->GetXaxis()->SetTitle(title);
0318   graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
0319   graph->Draw();
0320 
0321   //  graph->SetLineColor(kGray);
0322   //  graph->SetMarkerColor(kBlack);
0323   //  graph->Draw("AP");
0324 
0325   //   if( debug_ ) {
0326   //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
0327   //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
0328   //   }
0329   //   else {
0330   //     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
0331   //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
0332   //     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
0333   //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
0334   //   }
0335   //   sigmaPtVsEtaPlusErrTH1D->Draw();
0336   //   sigmaPtVsEtaTH1D->Draw("SAME");
0337   //   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
0338 
0339   sigmaPtVsEtaPlusErrTH1D->Write();
0340   sigmaPtVsEtaTH1D->Write();
0341   sigmaPtVsEtaMinusErrTH1D->Write();
0342 
0343   sigmaPtVsEtaPlusErr_->Write();
0344   sigmaPtVsEta_->Write();
0345   sigmaPtVsEtaMinusErr_->Write();
0346 
0347   // Mass
0348   sigmaMassVsEtaPlusErrTH1D->Write();
0349   sigmaMassVsEtaTH1D->Write();
0350   sigmaMassVsEtaMinusErrTH1D->Write();
0351 
0352   sigmaMassOverMassVsEtaPlusErrTH1D->Write();
0353   sigmaMassOverMassVsEtaTH1D->Write();
0354   sigmaMassOverMassVsEtaMinusErrTH1D->Write();
0355 
0356   sigmaMassVsEtaPlusErr_->Write();
0357   sigmaMassVsEta_->Write();
0358   sigmaMassVsEtaMinusErr_->Write();
0359 
0360   sigmaMassOverMassVsEtaPlusErr_->Write();
0361   sigmaMassOverMassVsEta_->Write();
0362   sigmaMassOverMassVsEtaMinusErr_->Write();
0363 
0364   canvas->Write();
0365 }
0366 
0367 // ------------ method called to for each event  ------------
0368 void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {}
0369 
0370 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
0371                                                  const lorentzVector& mu2,
0372                                                  const std::vector<double>& parval,
0373                                                  const double& sigmaPt1,
0374                                                  const double& sigmaPt2) {
0375   double* p = new double[(int)(parval.size())];
0376   std::vector<double>::const_iterator it = parval.begin();
0377   int id = 0;
0378   for (; it != parval.end(); ++it, ++id) {
0379     p[id] = *it;
0380   }
0381   double massRes = massResolution(mu1, mu2, p, sigmaPt1, sigmaPt2);
0382   delete[] p;
0383   return massRes;
0384 }
0385 
0386 double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
0387                                                  const lorentzVector& mu2,
0388                                                  double* parval,
0389                                                  const double& sigmaPt1,
0390                                                  const double& sigmaPt2) {
0391   double mass = (mu1 + mu2).mass();
0392   double pt1 = mu1.Pt();
0393   double phi1 = mu1.Phi();
0394   double eta1 = mu1.Eta();
0395   double theta1 = 2 * atan(exp(-eta1));
0396   double pt2 = mu2.Pt();
0397   double phi2 = mu2.Phi();
0398   double eta2 = mu2.Eta();
0399   double theta2 = 2 * atan(exp(-eta2));
0400 
0401   // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
0402   // -----------------------------------------------------------------------------------------------------------
0403   double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
0404                        sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
0405                             (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
0406                    pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
0407                   mass;
0408   double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
0409                        sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
0410                             (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
0411                    pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
0412                   mass;
0413   double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
0414   double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
0415   double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
0416                            sqrt((std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2) /
0417                                 (std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2)) -
0418                        pt1 * pt2 * cos(theta2) / sin(theta2)) /
0419                       mass;
0420   double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
0421                            sqrt((std::pow(pt1 / sin(theta1), 2) + MuScleFitUtils::mMu2) /
0422                                 (std::pow(pt2 / sin(theta2), 2) + MuScleFitUtils::mMu2)) -
0423                        pt2 * pt1 * cos(theta1) / sin(theta1)) /
0424                       mass;
0425 
0426   // Resolution parameters:
0427   // ----------------------
0428   // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
0429   // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
0430   double sigma_pt1 = sigmaPt1;
0431   double sigma_pt2 = sigmaPt2;
0432   double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt1, eta1, parval);
0433   double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi(pt2, eta2, parval);
0434   double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt1, eta1, parval);
0435   double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh(pt2, eta2, parval);
0436   double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2(pt1, eta1, pt2, eta2, parval);
0437 
0438   // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
0439   // multiply it by pt.
0440   double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
0441                          std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
0442                          std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2) +
0443                          2 * dmdpt1 * dmdpt2 * cov_pt1pt2);
0444 
0445   return mass_res;
0446 }
0447 
0448 void ErrorsPropagationAnalyzer::fillHistograms() {
0449   std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
0450   RootTreeHandler rootTreeHandler;
0451 
0452   typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
0453   MuonPairVector savedPair;
0454   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0455   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
0456   // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
0457 
0458   resolutionFunctionBase<std::vector<double> >* resolutionFunctionForVec = resolutionFunctionVecService(resolFitType_);
0459   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType_);
0460   MuScleFitUtils::debugMassResol_ = false;
0461   MuScleFitUtils::debug = 0;
0462   MuScleFitUtils::resfind = std::vector<int>(6, 0);
0463 
0464   SigmaPt sigmaPt(parameters_, errors_);
0465   SigmaPtDiff sigmaPtDiff;
0466 
0467   // Loop on all the pairs
0468   unsigned int i = 0;
0469   MuonPairVector::iterator it = savedPair.begin();
0470   std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
0471   for (; it != savedPair.end(); ++it, ++i) {
0472     double pt1 = it->first.pt();
0473     double eta1 = it->first.eta();
0474     double pt2 = it->second.pt();
0475     double eta2 = it->second.eta();
0476 
0477     if (debug_) {
0478       std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
0479     }
0480     // double fabsEta1 = fabs(eta1);
0481     // double fabsEta2 = fabs(eta2);
0482 
0483     if (pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0)
0484       continue;
0485 
0486     // double sigmaPt1 = sigmaPt( eta1 );
0487     // double sigmaPt2 = sigmaPt( eta2 );
0488 
0489     // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
0490     // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
0491     // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
0492     // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);
0493 
0494     double sigmaPt1 = resolutionFunctionForVec->sigmaPt(pt1, eta1, parameters_);
0495     double sigmaPt2 = resolutionFunctionForVec->sigmaPt(pt2, eta2, parameters_);
0496     double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
0497     double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
0498     double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError(pt1, eta1, parameters_, errors_);
0499     double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError(pt2, eta2, parameters_, errors_);
0500 
0501     // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
0502     // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
0503     // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
0504     double sigmaMass = massResolution(it->first, it->second, parameters_, sigmaPt1, sigmaPt2);
0505     double sigmaMassPlusErr = massResolution(it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2);
0506     double sigmaMassMinusErr = massResolution(it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2);
0507 
0508     double mass = (it->first + it->second).mass();
0509 
0510     if (debug_) {
0511       std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
0512       std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
0513       std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
0514     }
0515 
0516     // Protections from nans
0517     if (pt1 != pt1)
0518       continue;
0519     if (pt2 != pt2)
0520       continue;
0521     if (eta1 != eta1)
0522       continue;
0523     if (eta2 != eta2)
0524       continue;
0525     if (sigmaPt1 != sigmaPt1)
0526       continue;
0527     if (sigmaPt2 != sigmaPt2)
0528       continue;
0529     if (sigmaPtPlusErr1 != sigmaPtPlusErr1)
0530       continue;
0531     if (sigmaPtPlusErr2 != sigmaPtPlusErr2)
0532       continue;
0533     if (sigmaPtMinusErr1 != sigmaPtMinusErr1)
0534       continue;
0535     if (sigmaPtMinusErr2 != sigmaPtMinusErr2)
0536       continue;
0537     if (sigmaMass != sigmaMass)
0538       continue;
0539     if (sigmaMassPlusErr != sigmaMassPlusErr)
0540       continue;
0541     if (sigmaMassMinusErr != sigmaMassMinusErr)
0542       continue;
0543     if (mass == 0)
0544       continue;
0545 
0546     std::cout << "Muon pair number " << i << std::endl;
0547 
0548     // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
0549     // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
0550     // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
0551 
0552     if ((pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_)) {
0553       sigmaPtVsPt_->Fill(pt1, sigmaPt1);
0554       sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
0555       sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
0556       sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
0557       sigmaMassVsPt_->Fill(pt1, sigmaMass * mass);
0558       sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr * mass);
0559       sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr * mass);
0560       sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
0561       sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
0562       sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
0563 
0564       sigmaPtVsEta_->Fill(eta1, sigmaPt1);
0565       sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
0566       sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
0567       sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
0568       sigmaMassVsEta_->Fill(eta1, sigmaMass * mass);
0569       sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr * mass);
0570       sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr * mass);
0571       sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
0572       sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
0573       sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
0574     }
0575     if ((pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_)) {
0576       sigmaPtVsPt_->Fill(pt2, sigmaPt2);
0577       sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
0578       sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
0579       sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
0580       sigmaMassVsPt_->Fill(pt2, sigmaMass * mass);
0581       sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr * mass);
0582       sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr * mass);
0583       sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
0584       sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
0585       sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
0586 
0587       sigmaPtVsEta_->Fill(eta2, sigmaPt2);
0588       sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
0589       sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
0590       sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
0591       sigmaMassVsEta_->Fill(eta2, sigmaMass * mass);
0592       sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr * mass);
0593       sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr * mass);
0594       sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
0595       sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
0596       sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
0597     }
0598   }
0599 }
0600 
0601 //define this as a plug-in
0602 DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!