ErrorsPropagationAnalyzer.cc

Go to the documentation of this file.

#ifndef ERRORSANALYZER_CC
#define ERRORSANALYZER_CC

#include "ErrorsPropagationAnalyzer.h"

ErrorsPropagationAnalyzer::ErrorsPropagationAnalyzer(const edm::ParameterSet& iConfig) :
  treeFileName_( iConfig.getParameter<std::string>("InputFileName") ),
  resolFitType_( iConfig.getParameter<int>("ResolFitType") ),
  maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
  outputFileName_( iConfig.getParameter<std::string>("OutputFileName") ),
  ptBins_( iConfig.getParameter<int>("PtBins") ),
  ptMin_( iConfig.getParameter<double>("PtMin") ),
  ptMax_( iConfig.getParameter<double>("PtMax") ),
  etaBins_( iConfig.getParameter<int>("EtaBins") ),
  etaMin_( iConfig.getParameter<double>("EtaMin") ),
  etaMax_( iConfig.getParameter<double>("EtaMax") ),
  debug_( iConfig.getParameter<bool>("Debug") ),
  ptMinCut_( iConfig.getUntrackedParameter<double>("PtMinCut", 0.) ),
  ptMaxCut_( iConfig.getUntrackedParameter<double>("PtMaxCut", 999999.) ),
  etaMinCut_( iConfig.getUntrackedParameter<double>("EtaMinCut", 0.) ),
  etaMaxCut_( iConfig.getUntrackedParameter<double>("EtaMaxCut", 100.) )
{
  parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
  errors_ = iConfig.getParameter<std::vector<double> >("Errors");
  errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");

  if( (parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size()) ) {
    std::cout << "Error: parameters and errors have different number of values" << std::endl;
    exit(1);
  }

  fillValueError();

  sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
  sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
  sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);

  sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
  sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
  sigmaPtVsEtaMinusErr_ = new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);

  sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
  sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
  sigmaMassVsPtMinusErr_ = new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);

  sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
  sigmaMassVsEtaPlusErr_ = new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
  sigmaMassVsEtaMinusErr_ = new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);

  sigmaMassOverMassVsPt_ = new TProfile("sigmaMassOverMassVsPtProfile", "sigmaMassOverMassVsPt", ptBins_, ptMin_, ptMax_);
  sigmaMassOverMassVsPtPlusErr_ = new TProfile("sigmaMassOverMassVsPtPlusErrProfile", "sigmaMassOverMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
  sigmaMassOverMassVsPtMinusErr_ = new TProfile("sigmaMassOverMassVsPtMinusErrProfile", "sigmaMassOverMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);

  sigmaMassOverMassVsEta_ = new TProfile("sigmaMassOverMassVsEtaProfile", "sigmaMassOverMassVsEta", etaBins_, etaMin_, etaMax_);
  sigmaMassOverMassVsEtaPlusErr_ = new TProfile("sigmaMassOverMassVsEtaPlusErrProfile", "sigmaMassOverMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
  sigmaMassOverMassVsEtaMinusErr_ = new TProfile("sigmaMassOverMassVsEtaMinusErrProfile", "sigmaMassOverMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);

  sigmaPtVsPtDiff_ = new TProfile("sigmaPtVsPtDiffProfile", "sigmaPtVsPtDiff", ptBins_, ptMin_, ptMax_);
  sigmaPtVsEtaDiff_ = new TProfile("sigmaPtVsEtaDiffProfile", "sigmaPtVsEtaDiff", etaBins_, etaMin_, etaMax_);
}

void ErrorsPropagationAnalyzer::fillValueError()
{
  valuePlusError_.resize(parameters_.size());
  valueMinusError_.resize(parameters_.size());

  std::vector<double>::const_iterator parIt = parameters_.begin();
  std::vector<double>::const_iterator errIt = errors_.begin();
  std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
  int i=0;
  for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
    valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
    valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
  }
}

ErrorsPropagationAnalyzer::~ErrorsPropagationAnalyzer()
{
  gROOT->SetStyle("Plain");

  fillHistograms();

  TFile * outputFile = new TFile(outputFileName_, "RECREATE");

  drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta", "#sigma(Pt)/Pt");
  drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt", "#sigma(Pt)/Pt");

  drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta", "#sigma(M)");
  drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt", "#sigma(M)");

  drawHistograms(sigmaMassOverMassVsEta_, sigmaMassOverMassVsEtaPlusErr_, sigmaMassOverMassVsEtaMinusErr_, "sigmaMassOverMassVsEta", "#sigma(M)/M");
  drawHistograms(sigmaMassOverMassVsPt_, sigmaMassOverMassVsPtPlusErr_, sigmaMassOverMassVsPtMinusErr_, "sigmaMassOverMassVsPt", "#sigma(M)/M");

  sigmaPtVsPtDiff_->Write();
  sigmaPtVsEtaDiff_->Write();

  outputFile->Write();
  outputFile->Close();
}

void ErrorsPropagationAnalyzer::drawHistograms(const TProfile * histo, const TProfile * histoPlusErr,
                           const TProfile * histoMinusErr, const TString & type,
                           const TString & yLabel)
{
  TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
                     histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
                        histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
                         histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());

  TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
                       histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
                          histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
                           histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());

  TH1D * sigmaMassOverMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
                           histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaMassOverMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
                                  histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
  TH1D * sigmaMassOverMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
                                   histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());

  TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
  for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
    sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
    sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
    sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
  }
  int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
  // Draw TGraph with asymmetric errors
  Double_t * values = sigmaPtVsEtaTH1D->GetArray();
  Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
  Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
  double * posErrors = new double[numBins];
  double * negErrors = new double[numBins];

  TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
  TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);

  for( int i=1; i<=numBins; ++i ) {
    // std::cout << "filling " << i << std::endl;
    posErrors[i-1] = valuesPlus[i] - values[i];
    if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
    else negErrors[i-1] = values[i] - valuesMinus[i];

    graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
    graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
  }

  canvas->Draw();

  graphAsymmErrors->SetTitle("");
  graphAsymmErrors->SetFillColor(kGray);
  graphAsymmErrors->Draw("A2");
  TString title(type);
  if( type == "Eta" ) title = "#eta";
  graphAsymmErrors->GetXaxis()->SetTitle(title);
  graphAsymmErrors->GetYaxis()->SetTitle("yLabel");
  graph->Draw();

  //  graph->SetLineColor(kGray);
  //  graph->SetMarkerColor(kBlack);
  //  graph->Draw("AP");

//   if( debug_ ) {
//     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
//     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
//   }
//   else {
//     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
//     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
//     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
//     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
//   }
//   sigmaPtVsEtaPlusErrTH1D->Draw();
//   sigmaPtVsEtaTH1D->Draw("SAME");
//   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");

  sigmaPtVsEtaPlusErrTH1D->Write();
  sigmaPtVsEtaTH1D->Write();
  sigmaPtVsEtaMinusErrTH1D->Write();

  sigmaPtVsEtaPlusErr_->Write();
  sigmaPtVsEta_->Write();
  sigmaPtVsEtaMinusErr_->Write();

  // Mass
  sigmaMassVsEtaPlusErrTH1D->Write();
  sigmaMassVsEtaTH1D->Write();
  sigmaMassVsEtaMinusErrTH1D->Write();

  sigmaMassOverMassVsEtaPlusErrTH1D->Write();
  sigmaMassOverMassVsEtaTH1D->Write();
  sigmaMassOverMassVsEtaMinusErrTH1D->Write();

  sigmaMassVsEtaPlusErr_->Write();
  sigmaMassVsEta_->Write();
  sigmaMassVsEtaMinusErr_->Write();

  sigmaMassOverMassVsEtaPlusErr_->Write();
  sigmaMassOverMassVsEta_->Write();
  sigmaMassOverMassVsEtaMinusErr_->Write();

  canvas->Write();
}

// ------------ method called to for each event  ------------
void ErrorsPropagationAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
}

double ErrorsPropagationAnalyzer::massResolution(const lorentzVector& mu1,
                         const lorentzVector& mu2,
                         const std::vector< double >& parval,
                         const double& sigmaPt1,
                         const double& sigmaPt2)
{
  double * p = new double[(int)(parval.size())];
  std::vector<double>::const_iterator it = parval.begin();
  int id = 0;
  for ( ; it!=parval.end(); ++it, ++id) {
    p[id] = *it;
  }
  double massRes = massResolution (mu1, mu2, p, sigmaPt1, sigmaPt2);
  delete[] p;
  return massRes;
}

double ErrorsPropagationAnalyzer::massResolution( const lorentzVector& mu1,
                          const lorentzVector& mu2,
                          double* parval,
                          const double & sigmaPt1,
                          const double & sigmaPt2)
{
  double mass   = (mu1+mu2).mass();
  double pt1    = mu1.Pt();
  double phi1   = mu1.Phi();
  double eta1   = mu1.Eta();
  double theta1 = 2*atan(exp(-eta1));
  double pt2    = mu2.Pt();
  double phi2   = mu2.Phi();
  double eta2   = mu2.Eta();
  double theta2 = 2*atan(exp(-eta2));

  // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
  // -----------------------------------------------------------------------------------------------------------
  double dmdpt1  = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)/(std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2))-
            pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
  double dmdpt2  = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)/(std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2))-
            pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
  double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
  double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
  double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
                       sqrt((std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)/(std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)) -
               pt1*pt2*cos(theta2)/sin(theta2))/mass;
  double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
                       sqrt((std::pow(pt1/sin(theta1),2)+MuScleFitUtils::mMu2)/(std::pow(pt2/sin(theta2),2)+MuScleFitUtils::mMu2)) -
               pt2*pt1*cos(theta1)/sin(theta1))/mass;

  // Resolution parameters:
  // ----------------------
  // double sigma_pt1 = MuScleFitUtils::resolutionFunction->sigmaPt( pt1,eta1,parval );
  // double sigma_pt2 = MuScleFitUtils::resolutionFunction->sigmaPt( pt2,eta2,parval );
  double sigma_pt1 = sigmaPt1;
  double sigma_pt2 = sigmaPt2;
  double sigma_phi1 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt1,eta1,parval );
  double sigma_phi2 = MuScleFitUtils::resolutionFunction->sigmaPhi( pt2,eta2,parval );
  double sigma_cotgth1 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt1,eta1,parval );
  double sigma_cotgth2 = MuScleFitUtils::resolutionFunction->sigmaCotgTh( pt2,eta2,parval );
  double cov_pt1pt2 = MuScleFitUtils::resolutionFunction->covPt1Pt2( pt1, eta1, pt2, eta2, parval );

  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
  // multiply it by pt.
  double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2)+
             std::pow(dmdphi1*sigma_phi1,2)+std::pow(dmdphi2*sigma_phi2,2)+
             std::pow(dmdcotgth1*sigma_cotgth1,2)+std::pow(dmdcotgth2*sigma_cotgth2,2)+
             2*dmdpt1*dmdpt2*cov_pt1pt2);

  return mass_res;
}


void ErrorsPropagationAnalyzer::fillHistograms()
{
  std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
  RootTreeHandler rootTreeHandler;

  typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
  MuonPairVector savedPair;
  std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
  rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
  // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));

  resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec = resolutionFunctionVecService( resolFitType_ );
  MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType_ );
  MuScleFitUtils::debugMassResol_ = false;
  MuScleFitUtils::debug = 0;
  MuScleFitUtils::resfind = std::vector<int>(6, 0);

  SigmaPt sigmaPt( parameters_, errors_ );
  SigmaPtDiff sigmaPtDiff;

  // Loop on all the pairs
  unsigned int i = 0;
  MuonPairVector::iterator it = savedPair.begin();
  std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
  for( ; it != savedPair.end(); ++it, ++i ) {
    double pt1 = it->first.pt();
    double eta1 = it->first.eta();
    double pt2 = it->second.pt();
    double eta2 = it->second.eta();

    if( debug_ ) {
      std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
    }
    // double fabsEta1 = fabs(eta1);
    // double fabsEta2 = fabs(eta2);

    if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;


    // double sigmaPt1 = sigmaPt( eta1 );
    // double sigmaPt2 = sigmaPt( eta2 );

    // double sigmaPtPlusErr1 = sigmaPt1 + sigmaPt.sigma(eta1);
    // double sigmaPtPlusErr2 = sigmaPt2 + sigmaPt.sigma(eta2);
    // double sigmaPtMinusErr1 = sigmaPt1 - sigmaPt.sigma(eta1);
    // double sigmaPtMinusErr2 = sigmaPt2 - sigmaPt.sigma(eta2);

    double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
    double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
    double sigmaPtPlusErr1 = sigmaPt1 + resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_);
    double sigmaPtPlusErr2 = sigmaPt2 + resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );
    double sigmaPtMinusErr1 = sigmaPt1 - resolutionFunctionForVec->sigmaPtError( pt1,eta1,parameters_,errors_ );
    double sigmaPtMinusErr2 = sigmaPt2 - resolutionFunctionForVec->sigmaPtError( pt2,eta2,parameters_,errors_ );

    // double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
    // double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
    // double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
    double sigmaMass = massResolution( it->first, it->second, parameters_, sigmaPt1, sigmaPt2 );
    double sigmaMassPlusErr = massResolution( it->first, it->second, parameters_, sigmaPtPlusErr1, sigmaPtPlusErr2 );
    double sigmaMassMinusErr = massResolution( it->first, it->second, parameters_, sigmaPtMinusErr1, sigmaPtMinusErr2 );

    double mass = (it->first + it->second).mass();

    if( debug_ ) {
      std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
      std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
      std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
    }

    // Protections from nans
    if( pt1 != pt1 ) continue;
    if( pt2 != pt2 ) continue;
    if( eta1 != eta1 ) continue;
    if( eta2 != eta2 ) continue;
    if( sigmaPt1 != sigmaPt1 ) continue;
    if( sigmaPt2 != sigmaPt2 ) continue;
    if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
    if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
    if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
    if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
    if( sigmaMass != sigmaMass ) continue;
    if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
    if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
    if( mass == 0 ) continue;

    std::cout << "Muon pair number " << i << std::endl;

    // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
    // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
    // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;


    if( (pt1 >= ptMinCut_ && pt1 <= ptMaxCut_) && (fabs(eta1) >= etaMinCut_ && fabs(eta1) <= etaMaxCut_) ) {
      sigmaPtVsPt_->Fill(pt1, sigmaPt1);
      sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
      sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
      sigmaPtVsPtDiff_->Fill(pt1, sigmaPtDiff.squaredDiff(eta1));
      sigmaMassVsPt_->Fill(pt1, sigmaMass*mass);
      sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr*mass);
      sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr*mass);
      sigmaMassOverMassVsPt_->Fill(pt1, sigmaMass);
      sigmaMassOverMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
      sigmaMassOverMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);

      sigmaPtVsEta_->Fill(eta1, sigmaPt1);
      sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
      sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
      sigmaPtVsEtaDiff_->Fill(eta1, sigmaPtDiff.squaredDiff(eta1));
      sigmaMassVsEta_->Fill(eta1, sigmaMass*mass);
      sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr*mass);
      sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr*mass);
      sigmaMassOverMassVsEta_->Fill(eta1, sigmaMass);
      sigmaMassOverMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
      sigmaMassOverMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
    }
    if( (pt2 >= ptMinCut_ && pt2 <= ptMaxCut_) && (fabs(eta2) >= etaMinCut_ && fabs(eta2) <= etaMaxCut_) ) {
      sigmaPtVsPt_->Fill(pt2, sigmaPt2);
      sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
      sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
      sigmaPtVsPtDiff_->Fill(pt2, sigmaPtDiff.squaredDiff(eta2));
      sigmaMassVsPt_->Fill(pt2, sigmaMass*mass);
      sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr*mass);
      sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr*mass);
      sigmaMassOverMassVsPt_->Fill(pt2, sigmaMass);
      sigmaMassOverMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
      sigmaMassOverMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);

      sigmaPtVsEta_->Fill(eta2, sigmaPt2);
      sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
      sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
      sigmaPtVsEtaDiff_->Fill(eta2, sigmaPtDiff.squaredDiff(eta2));
      sigmaMassVsEta_->Fill(eta2, sigmaMass*mass);
      sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr*mass);
      sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr*mass);
      sigmaMassOverMassVsEta_->Fill(eta2, sigmaMass);
      sigmaMassOverMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
      sigmaMassOverMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
    }
  }
}

//define this as a plug-in
DEFINE_FWK_MODULE(ErrorsPropagationAnalyzer);

#endif