Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:41 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ResolutionAnalyzer
0004 // Class:      ResolutionAnalyzer
0005 //
0006 /**\class ResolutionAnalyzer ResolutionAnalyzer.cc MuonAnalysis/MomentumScaleCalibration/plugins/ResolutionAnalyzer.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 #include <CLHEP/Vector/LorentzVector.h>
0024 
0025 // user include files
0026 #include "DataFormats/Candidate/interface/Candidate.h"
0027 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0028 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0029 #include "DataFormats/MuonReco/interface/Muon.h"
0030 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "HepMC/GenEvent.h"
0038 #include "HepMC/GenParticle.h"
0039 #include "HepPDT/ParticleDataTable.hh"
0040 #include "HepPDT/TableBuilder.hh"
0041 #include "HepPDT/defs.h"
0042 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
0043 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
0044 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0045 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0046 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0047 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0048 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0049 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0050 
0051 #include "Histograms.h"
0052 #include "MuScleFitUtils.h"
0053 //
0054 // class decleration
0055 //
0056 
0057 class ResolutionAnalyzer : public edm::one::EDAnalyzer<> {
0058 public:
0059   explicit ResolutionAnalyzer(const edm::ParameterSet&);
0060   ~ResolutionAnalyzer() override;
0061 
0062 private:
0063   void analyze(const edm::Event&, const edm::EventSetup&) override;
0064   void endJob() override{};
0065 
0066   template <typename T>
0067   std::vector<reco::LeafCandidate> fillMuonCollection(const std::vector<T>& tracks) {
0068     std::vector<reco::LeafCandidate> muons;
0069     typename std::vector<T>::const_iterator track;
0070     for (track = tracks.begin(); track != tracks.end(); ++track) {
0071       reco::Particle::LorentzVector mu(
0072           track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
0073       MuScleFitUtils::goodmuon++;
0074       if (debug_ > 0)
0075         std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon
0076                   << ": initial value   Pt = " << mu.Pt() << std::endl;
0077       reco::LeafCandidate muon(track->charge(), mu);
0078       // Store muon
0079       // ----------
0080       muons.push_back(muon);
0081     }
0082     return muons;
0083   }
0084 
0085   /// Used to fill the map with the histograms needed
0086   void fillHistoMap();
0087   /// Writes the histograms in the map
0088   void writeHistoMap();
0089   /// Returns true if the two particles have DeltaR < cut
0090   bool checkDeltaR(const reco::Particle::LorentzVector& genMu, const reco::Particle::LorentzVector& recMu);
0091 
0092   // ----------member data ---------------------------
0093 
0094   // Collections labels
0095   // ------------------
0096   edm::InputTag theMuonLabel_;
0097 
0098   int theMuonType_;
0099   std::string theRootFileName_;
0100   std::string theCovariancesRootFileName_;
0101   bool debug_;
0102   std::map<std::string, Histograms*> mapHisto_;
0103   TFile* outputFile_;
0104 
0105   int eventCounter_;
0106   bool resonance_;
0107   bool readCovariances_;
0108 
0109   TString treeFileName_;
0110   int32_t maxEvents_;
0111 
0112   double ptMax_;
0113 
0114   HCovarianceVSxy* massResolutionVsPtEta_;
0115   TH2D* recoPtVsgenPt_;
0116   TH2D* recoPtVsgenPtEta12_;
0117   TH1D* deltaPtOverPt_;
0118   TH1D* deltaPtOverPtForEta12_;
0119 };
0120 
0121 //
0122 // constructors and destructor
0123 //
0124 ResolutionAnalyzer::ResolutionAnalyzer(const edm::ParameterSet& iConfig)
0125     : theMuonLabel_(iConfig.getParameter<edm::InputTag>("MuonLabel")),
0126       theMuonType_(iConfig.getParameter<int>("MuonType")),
0127       theRootFileName_(iConfig.getUntrackedParameter<std::string>("OutputFileName")),
0128       theCovariancesRootFileName_(iConfig.getUntrackedParameter<std::string>("InputFileName")),
0129       debug_(iConfig.getUntrackedParameter<bool>("Debug")),
0130       resonance_(iConfig.getParameter<bool>("Resonance")),
0131       readCovariances_(iConfig.getParameter<bool>("ReadCovariances")),
0132       treeFileName_(iConfig.getParameter<std::string>("InputTreeName")),
0133       maxEvents_(iConfig.getParameter<int>("MaxEvents")),
0134       ptMax_(iConfig.getParameter<double>("PtMax")) {
0135   //now do what ever initialization is needed
0136 
0137   // Initial parameters values
0138   // -------------------------
0139   int resolFitType = iConfig.getParameter<int>("ResolFitType");
0140   MuScleFitUtils::ResolFitType = resolFitType;
0141   // MuScleFitUtils::resolutionFunction = resolutionFunctionArray[resolFitType];
0142   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType);
0143   // MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionArrayForVec[resolFitType];
0144   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolFitType);
0145 
0146   MuScleFitUtils::parResol = iConfig.getParameter<std::vector<double> >("parResol");
0147 
0148   MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("ResFind");
0149 
0150   outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
0151   outputFile_->cd();
0152   fillHistoMap();
0153 
0154   eventCounter_ = 0;
0155 }
0156 
0157 ResolutionAnalyzer::~ResolutionAnalyzer() {
0158   outputFile_->cd();
0159   writeHistoMap();
0160   outputFile_->Close();
0161   std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
0162 }
0163 
0164 //
0165 // member functions
0166 //
0167 
0168 // ------------ method called to for each event  ------------
0169 void ResolutionAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0170   using namespace edm;
0171 
0172   std::cout << "starting" << std::endl;
0173 
0174   lorentzVector nullLorentzVector(0, 0, 0, 0);
0175 
0176   RootTreeHandler rootTreeHandler;
0177   typedef std::vector<std::pair<lorentzVector, lorentzVector> > MuonPairVector;
0178   MuonPairVector savedPairVector;
0179   MuonPairVector genPairVector;
0180 
0181   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0182   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPairVector, 0, &evtRun, &genPairVector);
0183   MuonPairVector::iterator savedPair = savedPairVector.begin();
0184   MuonPairVector::iterator genPair = genPairVector.begin();
0185   std::cout << "Starting loop on " << savedPairVector.size() << " muons" << std::endl;
0186   for (; savedPair != savedPairVector.end(); ++savedPair, ++genPair) {
0187     ++eventCounter_;
0188 
0189     if ((eventCounter_ % 10000) == 0) {
0190       std::cout << "event = " << eventCounter_ << std::endl;
0191     }
0192 
0193     lorentzVector recMu1(savedPair->first);
0194     lorentzVector recMu2(savedPair->second);
0195 
0196     if (resonance_) {
0197       // Histograms with genParticles characteristics
0198       // --------------------------------------------
0199 
0200       reco::Particle::LorentzVector genMother(genPair->first + genPair->second);
0201 
0202       mapHisto_["GenMother"]->Fill(genMother);
0203       mapHisto_["DeltaGenMotherMuons"]->Fill(genPair->first, genPair->second);
0204       mapHisto_["GenMotherMuons"]->Fill(genPair->first);
0205       mapHisto_["GenMotherMuons"]->Fill(genPair->second);
0206 
0207       // Match the reco muons with the gen and sim tracks
0208       // ------------------------------------------------
0209       if (checkDeltaR(genPair->first, recMu1)) {
0210         mapHisto_["PtResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Pt() + recMu1.Pt()) / genPair->first.Pt(), -1);
0211         mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Theta() + recMu1.Theta()), -1);
0212         mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(
0213             recMu1,
0214             (-cos(genPair->first.Theta()) / sin(genPair->first.Theta()) + cos(recMu1.Theta()) / sin(recMu1.Theta())),
0215             -1);
0216         mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu1, (-genPair->first.Eta() + recMu1.Eta()), -1);
0217         // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu1,(-genPair->first.Phi()+recMu1.Phi()),-1);
0218         mapHisto_["PhiResolutionGenVSMu"]->Fill(
0219             recMu1, MuScleFitUtils::deltaPhiNoFabs(recMu1.Phi(), genPair->first.Phi()), -1);
0220         recoPtVsgenPt_->Fill(genPair->first.Pt(), recMu1.Pt());
0221         deltaPtOverPt_->Fill((recMu1.Pt() - genPair->first.Pt()) / genPair->first.Pt());
0222         if (fabs(recMu1.Eta()) > 1 && fabs(recMu1.Eta()) < 1.2) {
0223           recoPtVsgenPtEta12_->Fill(genPair->first.Pt(), recMu1.Pt());
0224           deltaPtOverPtForEta12_->Fill((recMu1.Pt() - genPair->first.Pt()) / genPair->first.Pt());
0225         }
0226       }
0227       if (checkDeltaR(genPair->second, recMu2)) {
0228         mapHisto_["PtResolutionGenVSMu"]->Fill(
0229             recMu2, (-genPair->second.Pt() + recMu2.Pt()) / genPair->second.Pt(), +1);
0230         mapHisto_["ThetaResolutionGenVSMu"]->Fill(recMu2, (-genPair->second.Theta() + recMu2.Theta()), +1);
0231         mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(
0232             recMu2,
0233             (-cos(genPair->second.Theta()) / sin(genPair->second.Theta()) + cos(recMu2.Theta()) / sin(recMu2.Theta())),
0234             +1);
0235         mapHisto_["EtaResolutionGenVSMu"]->Fill(recMu2, (-genPair->second.Eta() + recMu2.Eta()), +1);
0236         // mapHisto_["PhiResolutionGenVSMu"]->Fill(recMu2,(-genPair->second.Phi()+recMu2.Phi()),+1);
0237         mapHisto_["PhiResolutionGenVSMu"]->Fill(
0238             recMu2, MuScleFitUtils::deltaPhiNoFabs(recMu2.Phi(), genPair->second.Phi()), +1);
0239         recoPtVsgenPt_->Fill(genPair->second.Pt(), recMu2.Pt());
0240         deltaPtOverPt_->Fill((recMu2.Pt() - genPair->second.Pt()) / genPair->second.Pt());
0241         if (fabs(recMu2.Eta()) > 1 && fabs(recMu2.Eta()) < 1.2) {
0242           recoPtVsgenPtEta12_->Fill(genPair->second.Pt(), recMu2.Pt());
0243           deltaPtOverPtForEta12_->Fill((recMu2.Pt() - genPair->second.Pt()) / genPair->second.Pt());
0244         }
0245       }
0246 
0247       // Fill the mass resolution histograms
0248       // -----------------------------------
0249       // check if the recoMuons match the genMuons
0250       // if( MuScleFitUtils::ResFound && checkDeltaR(simMu.first,recMu1) && checkDeltaR(simMu.second,recMu2) ) {
0251       if (genPair->first != nullLorentzVector && genPair->second != nullLorentzVector &&
0252           checkDeltaR(genPair->first, recMu1) && checkDeltaR(genPair->second, recMu2)) {
0253         double recoMass = (recMu1 + recMu2).mass();
0254         double genMass = (genPair->first + genPair->second).mass();
0255         // first is always mu-, second is always mu+
0256         mapHisto_["MassResolution"]->Fill(recMu1, -1, genPair->first, recMu2, +1, genPair->second, recoMass, genMass);
0257 
0258         // Fill the reconstructed resonance
0259         reco::Particle::LorentzVector recoResonance(recMu1 + recMu2);
0260         mapHisto_["RecoResonance"]->Fill(recoResonance);
0261         mapHisto_["DeltaRecoResonanceMuons"]->Fill(recMu1, recMu2);
0262         mapHisto_["RecoResonanceMuons"]->Fill(recMu1);
0263         mapHisto_["RecoResonanceMuons"]->Fill(recMu2);
0264 
0265         // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
0266         if (genMass != 0) {
0267           // double diffMass = (recoMass - genMass)/genMass;
0268           double diffMass = recoMass - genMass;
0269           // Fill if for both muons
0270           double pt1 = recMu1.pt();
0271           double eta1 = recMu1.eta();
0272           double pt2 = recMu2.pt();
0273           double eta2 = recMu2.eta();
0274           // This is to avoid nan
0275           if (diffMass == diffMass) {
0276             massResolutionVsPtEta_->Fill(pt1, eta1, diffMass, diffMass);
0277             massResolutionVsPtEta_->Fill(pt2, eta2, diffMass, diffMass);
0278           } else {
0279             std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
0280           }
0281           // Fill with mass resolution from resolution function
0282           double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
0283           // The value given by massRes is already divided by the mass, since the derivative functions have mass at the denominator.
0284           // We alos take the squared value, since var = sigma^2.
0285           mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
0286           mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
0287         }
0288 
0289         // Fill resolution functions for the muons (fill the squared value to make it comparable with the variance)
0290         mapHisto_["hFunctionResolPt"]->Fill(
0291             recMu1,
0292             MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
0293             -1);
0294         mapHisto_["hFunctionResolCotgTheta"]->Fill(
0295             recMu1,
0296             MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
0297             -1);
0298         mapHisto_["hFunctionResolPhi"]->Fill(
0299             recMu1,
0300             MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), MuScleFitUtils::parResol),
0301             -1);
0302         mapHisto_["hFunctionResolPt"]->Fill(
0303             recMu2,
0304             MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
0305             +1);
0306         mapHisto_["hFunctionResolCotgTheta"]->Fill(
0307             recMu2,
0308             MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
0309             +1);
0310         mapHisto_["hFunctionResolPhi"]->Fill(
0311             recMu2,
0312             MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), MuScleFitUtils::parResol),
0313             +1);
0314 
0315         if (readCovariances_) {
0316           // Compute mass error terms
0317           // ------------------------
0318           double mass = (recMu1 + recMu2).mass();
0319           double pt1 = recMu1.Pt();
0320           double phi1 = recMu1.Phi();
0321           double eta1 = recMu1.Eta();
0322           double theta1 = 2 * atan(exp(-eta1));
0323           double pt2 = recMu2.Pt();
0324           double phi2 = recMu2.Phi();
0325           double eta2 = recMu2.Eta();
0326           double theta2 = 2 * atan(exp(-eta2));
0327           // Derivatives
0328           double mMu2 = MuScleFitUtils::mMu2;
0329           double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
0330                                sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0331                            pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
0332                           mass;
0333           double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
0334                                sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0335                            pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
0336                           mass;
0337           double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
0338           double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
0339           double dmdcotgth1 =
0340               (pt1 * pt1 * cos(theta1) / sin(theta1) *
0341                    sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0342                pt1 * pt2 * cos(theta2) / sin(theta2)) /
0343               mass;
0344           double dmdcotgth2 =
0345               (pt2 * pt2 * cos(theta2) / sin(theta2) *
0346                    sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0347                pt2 * pt1 * cos(theta1) / sin(theta1)) /
0348               mass;
0349 
0350           // Multiplied by the pt here
0351           // -------------------------
0352           double dmdpt[2] = {dmdpt1 * recMu1.Pt(), dmdpt2 * recMu2.Pt()};
0353           double dmdphi[2] = {dmdphi1, dmdphi2};
0354           double dmdcotgth[2] = {dmdcotgth1, dmdcotgth2};
0355 
0356           // Evaluate the single terms in the mass error expression
0357 
0358           reco::Particle::LorentzVector* recMu[2] = {&recMu1, &recMu2};
0359           int charge[2] = {-1, +1};
0360 
0361           double fullMassRes = 0.;
0362           double massRes1 = 0.;
0363           double massRes2 = 0.;
0364           double massRes3 = 0.;
0365           double massRes4 = 0.;
0366           double massRes5 = 0.;
0367           double massRes6 = 0.;
0368           double massResPtAndPt12 = 0.;
0369 
0370           for (int i = 0; i < 2; ++i) {
0371             double ptVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt");
0372             double cotgThetaVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta");
0373             double phiVariance = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi");
0374             double pt_cotgTheta = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-CotgTheta");
0375             double pt_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt-Phi");
0376             double cotgTheta_phi = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta-Phi");
0377 
0378             double pt1_pt2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt1-Pt2");
0379             double cotgTheta1_cotgTheta2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta1-CotgTheta2");
0380             double phi1_phi2 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Phi1-Phi2");
0381             double pt12_cotgTheta21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-CotgTheta21");
0382             double pt12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "Pt12-Phi21");
0383             double cotgTheta12_phi21 = mapHisto_["ReadCovariances"]->Get(*(recMu[i]), "CotgTheta12-Phi21");
0384 
0385             // ATTENTION: Pt covariance terms are multiplied by Pt, since DeltaPt/Pt was used to compute them
0386             mapHisto_["MassResolutionPt"]->Fill(*(recMu[i]), ptVariance * std::pow(dmdpt[i], 2), charge[i]);
0387             mapHisto_["MassResolutionCotgTheta"]->Fill(
0388                 *(recMu[i]), cotgThetaVariance * std::pow(dmdcotgth[i], 2), charge[i]);
0389             mapHisto_["MassResolutionPhi"]->Fill(*(recMu[i]), phiVariance * std::pow(dmdphi[i], 2), charge[i]);
0390             mapHisto_["MassResolutionPt-CotgTheta"]->Fill(
0391                 *(recMu[i]), 2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i], charge[i]);
0392             mapHisto_["MassResolutionPt-Phi"]->Fill(*(recMu[i]), 2 * pt_phi * dmdpt[i] * dmdphi[i], charge[i]);
0393             mapHisto_["MassResolutionCotgTheta-Phi"]->Fill(
0394                 *(recMu[i]), 2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i], charge[i]);
0395 
0396             mapHisto_["MassResolutionPt1-Pt2"]->Fill(*(recMu[i]), pt1_pt2 * dmdpt[0] * dmdpt[1], charge[i]);
0397             mapHisto_["MassResolutionCotgTheta1-CotgTheta2"]->Fill(
0398                 *(recMu[i]), cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1], charge[i]);
0399             mapHisto_["MassResolutionPhi1-Phi2"]->Fill(*(recMu[i]), phi1_phi2 * dmdphi[0] * dmdphi[1], charge[i]);
0400             // This must be filled for both configurations: 12 and 21
0401             mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill(
0402                 *(recMu[i]), pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1], charge[i]);
0403             mapHisto_["MassResolutionPt12-CotgTheta21"]->Fill(
0404                 *(recMu[i]), pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0], charge[i]);
0405             mapHisto_["MassResolutionPt12-Phi21"]->Fill(*(recMu[i]), pt12_phi21 * dmdpt[0] * dmdphi[1], charge[i]);
0406             mapHisto_["MassResolutionPt12-Phi21"]->Fill(*(recMu[i]), pt12_phi21 * dmdpt[1] * dmdphi[0], charge[i]);
0407             mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill(
0408                 *(recMu[i]), cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1], charge[i]);
0409             mapHisto_["MassResolutionCotgTheta12-Phi21"]->Fill(
0410                 *(recMu[i]), cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0], charge[i]);
0411 
0412             // Sigmas for comparison
0413             mapHisto_["sigmaPtFromVariance"]->Fill(*(recMu[i]), sqrt(ptVariance), charge[i]);
0414             mapHisto_["sigmaCotgThetaFromVariance"]->Fill(*(recMu[i]), sqrt(cotgThetaVariance), charge[i]);
0415             mapHisto_["sigmaPhiFromVariance"]->Fill(*(recMu[i]), sqrt(phiVariance), charge[i]);
0416 
0417             // Pt term from function
0418             mapHisto_["MassResolutionPtFromFunction"]->Fill(
0419                 *(recMu[i]),
0420                 (MuScleFitUtils::resolutionFunctionForVec->sigmaPt(
0421                     (recMu[i])->Pt(), (recMu[i])->Eta(), MuScleFitUtils::parResol)) *
0422                     std::pow(dmdpt[i], 2),
0423                 charge[i]);
0424 
0425             fullMassRes += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
0426                            phiVariance * std::pow(dmdphi[i], 2) +
0427 
0428                            // These are worth twice the others since there are: pt1-phi1, phi1-pt1, pt2-phi2, phi2-pt2
0429                            2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
0430                            2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
0431 
0432                            pt1_pt2 * dmdpt[0] * dmdpt[1] + cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] +
0433                            phi1_phi2 * dmdphi[0] * dmdphi[1] +
0434 
0435                            // These are filled twice, because of the two combinations
0436                            pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
0437                            pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
0438                            cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
0439 
0440             massRes1 += ptVariance * std::pow(dmdpt[i], 2);
0441             massRes2 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2);
0442             massRes3 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
0443                         phiVariance * std::pow(dmdphi[i], 2);
0444             massRes4 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
0445                         phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
0446                         2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
0447                         2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i];
0448             massRes5 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
0449                         phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
0450                         2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
0451                         2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
0452                         cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1];
0453             massRes6 += ptVariance * std::pow(dmdpt[i], 2) + cotgThetaVariance * std::pow(dmdcotgth[i], 2) +
0454                         phiVariance * std::pow(dmdphi[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1] +
0455                         2 * pt_cotgTheta * dmdpt[i] * dmdcotgth[i] + 2 * pt_phi * dmdpt[i] * dmdphi[i] +
0456                         2 * cotgTheta_phi * dmdcotgth[i] * dmdphi[i] +
0457                         cotgTheta1_cotgTheta2 * dmdcotgth[0] * dmdcotgth[1] + phi1_phi2 * dmdphi[0] * dmdphi[1] +
0458                         pt12_cotgTheta21 * dmdpt[0] * dmdcotgth[1] + pt12_cotgTheta21 * dmdpt[1] * dmdcotgth[0] +
0459                         pt12_phi21 * dmdpt[0] * dmdphi[1] + pt12_phi21 * dmdpt[1] * dmdphi[0] +
0460                         cotgTheta12_phi21 * dmdcotgth[0] * dmdphi[1] + cotgTheta12_phi21 * dmdcotgth[1] * dmdphi[0];
0461 
0462             massResPtAndPt12 += ptVariance * std::pow(dmdpt[i], 2) + pt1_pt2 * dmdpt[0] * dmdpt[1];
0463 
0464             // Derivatives
0465             mapHisto_["DerivativePt"]->Fill(*(recMu[i]), dmdpt[i], charge[i]);
0466             mapHisto_["DerivativeCotgTheta"]->Fill(*(recMu[i]), dmdcotgth[i], charge[i]);
0467             mapHisto_["DerivativePhi"]->Fill(*(recMu[i]), dmdphi[i], charge[i]);
0468           }
0469           // Fill the complete resolution function with covariance terms
0470           mapHisto_["FullMassResolution"]->Fill(*(recMu[0]), fullMassRes, charge[0]);
0471           mapHisto_["FullMassResolution"]->Fill(*(recMu[1]), fullMassRes, charge[1]);
0472 
0473           mapHisto_["MassRes1"]->Fill(*(recMu[0]), massRes1, charge[0]);
0474           mapHisto_["MassRes1"]->Fill(*(recMu[1]), massRes1, charge[1]);
0475           mapHisto_["MassRes2"]->Fill(*(recMu[0]), massRes2, charge[0]);
0476           mapHisto_["MassRes2"]->Fill(*(recMu[1]), massRes2, charge[1]);
0477           mapHisto_["MassRes3"]->Fill(*(recMu[0]), massRes3, charge[0]);
0478           mapHisto_["MassRes3"]->Fill(*(recMu[1]), massRes3, charge[1]);
0479           mapHisto_["MassRes4"]->Fill(*(recMu[0]), massRes4, charge[0]);
0480           mapHisto_["MassRes4"]->Fill(*(recMu[1]), massRes4, charge[1]);
0481           mapHisto_["MassRes5"]->Fill(*(recMu[0]), massRes5, charge[0]);
0482           mapHisto_["MassRes5"]->Fill(*(recMu[1]), massRes5, charge[1]);
0483           mapHisto_["MassRes6"]->Fill(*(recMu[0]), massRes6, charge[0]);
0484           mapHisto_["MassRes6"]->Fill(*(recMu[1]), massRes6, charge[1]);
0485           mapHisto_["MassResPtAndPt12"]->Fill(*(recMu[0]), massResPtAndPt12, charge[0]);
0486           mapHisto_["MassResPtAndPt12"]->Fill(*(recMu[1]), massResPtAndPt12, charge[1]);
0487         } else {
0488           // Fill the covariances histograms
0489           mapHisto_["Covariances"]->Fill(recMu1, genPair->first, recMu2, genPair->second);
0490         }
0491       }
0492     }  // end if resonance
0493   }
0494   //   else {
0495   //
0496   //     // Loop on the recMuons
0497   //     std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
0498   //     for ( ; recMuon!=muons.end(); ++recMuon ) {
0499   //       int charge = recMuon->charge();
0500   //
0501   //       lorentzVector recMu(recMuon->p4());
0502   //
0503   //       // Find the matching MC muon
0504   //       const HepMC::GenEvent* Evt = evtMC->GetEvent();
0505   //       //Loop on generated particles
0506   //       std::map<double, lorentzVector> genAssocMap;
0507   //       HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin();
0508   //       for( ; part!=Evt->particles_end(); ++part ) {
0509   //         if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
0510   //           lorentzVector genMu = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
0511   //                                                (*part)->momentum().pz(),(*part)->momentum().e()));
0512   //
0513   //           double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genPair->Phi()) +
0514   //                                ((recMu.Eta()-genPair->Eta()) * (recMu.Eta()-genPair->Eta())));
0515   //
0516   //           // 13 for the muon (-1) and -13 for the antimuon (+1), thus pdg*charge = -13.
0517   //           // Only in this case we consider it matching.
0518   //           if( ((*part)->pdg_id())*charge == -13 ) genAssocMap.insert(std::make_pair(deltaR, genMu));
0519   //         }
0520   //       }
0521   //       // Take the closest in deltaR
0522   //       lorentzVector genMu(genAssocMap.begin()->second);
0523   //
0524   //       // Histograms with genParticles characteristics
0525   //       // --------------------------------------------
0526   //
0527   //       if(checkDeltaR(genMu,recMu)){
0528   //         mapHisto_["PtResolutionGenVSMu"]->Fill(genMu,(-genPair->Pt()+recMu.Pt())/genPair->Pt(),charge);
0529   //         mapHisto_["ThetaResolutionGenVSMu"]->Fill(genMu,(-genPair->Theta()+recMu.Theta()),charge);
0530   //         mapHisto_["CotgThetaResolutionGenVSMu"]->Fill(genMu,(-cos(genPair->Theta())/sin(genPair->Theta())
0531   //                                                              +cos(recMu.Theta())/sin(recMu.Theta())),charge);
0532   //         mapHisto_["EtaResolutionGenVSMu"]->Fill(genMu,(-genPair->Eta()+recMu.Eta()),charge);
0533   //         mapHisto_["PhiResolutionGenVSMu"]->Fill(genMu,MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genPair->Phi()),charge);
0534   //       }
0535   //     }
0536 }
0537 
0538 void ResolutionAnalyzer::fillHistoMap() {
0539   outputFile_->cd();
0540 
0541   // Resonances
0542   // If no Z is required, use a smaller mass range.
0543   double minMass = 0.;
0544   double maxMass = 200.;
0545   if (MuScleFitUtils::resfind[0] != 1) {
0546     maxMass = 30.;
0547   }
0548   mapHisto_["GenMother"] = new HParticle(outputFile_, "GenMother", minMass, maxMass);
0549   mapHisto_["SimResonance"] = new HParticle(outputFile_, "SimResonance", minMass, maxMass);
0550   mapHisto_["RecoResonance"] = new HParticle(outputFile_, "RecoResonance", minMass, maxMass);
0551 
0552   // Resonance muons
0553   mapHisto_["GenMotherMuons"] = new HParticle(outputFile_, "GenMotherMuons", minMass, 1.);
0554   mapHisto_["SimResonanceMuons"] = new HParticle(outputFile_, "SimResonanceMuons", minMass, 1.);
0555   mapHisto_["RecoResonanceMuons"] = new HParticle(outputFile_, "RecoResonanceMuons", minMass, 1.);
0556 
0557   // Deltas between resonance muons
0558   mapHisto_["DeltaGenMotherMuons"] = new HDelta(outputFile_, "DeltaGenMotherMuons");
0559   mapHisto_["DeltaSimResonanceMuons"] = new HDelta(outputFile_, "DeltaSimResonanceMuons");
0560   mapHisto_["DeltaRecoResonanceMuons"] = new HDelta(outputFile_, "DeltaRecoResonanceMuons");
0561 
0562   //   //Reconstructed muon kinematics
0563   //   //-----------------------------
0564   //   mapHisto_["hRecBestMu"]             = new HParticle         ("hRecBestMu");
0565   //   mapHisto_["hRecBestMu_Acc"]         = new HParticle         ("hRecBestMu_Acc");
0566   //   mapHisto_["hDeltaRecBestMu"]        = new HDelta            ("hDeltaRecBestMu");
0567 
0568   //   mapHisto_["hRecBestRes"]            = new HParticle         ("hRecBestRes");
0569   //   mapHisto_["hRecBestRes_Acc"]        = new HParticle         ("hRecBestRes_Acc");
0570   //   mapHisto_["hRecBestResVSMu"]        = new HMassVSPart       ("hRecBestResVSMu");
0571 
0572   //Resolution VS muon kinematic
0573   //----------------------------
0574   mapHisto_["PtResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "PtResolutionGenVSMu");
0575   mapHisto_["PtResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "PtResolutionSimVSMu");
0576   mapHisto_["EtaResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "EtaResolutionGenVSMu");
0577   mapHisto_["EtaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "EtaResolutionSimVSMu");
0578   mapHisto_["ThetaResolutionGenVSMu"] = new HResolutionVSPart(outputFile_, "ThetaResolutionGenVSMu");
0579   mapHisto_["ThetaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "ThetaResolutionSimVSMu");
0580   mapHisto_["CotgThetaResolutionGenVSMu"] =
0581       new HResolutionVSPart(outputFile_, "CotgThetaResolutionGenVSMu", -0.02, 0.02, -0.02, 0.02);
0582   mapHisto_["CotgThetaResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "CotgThetaResolutionSimVSMu");
0583   mapHisto_["PhiResolutionGenVSMu"] =
0584       new HResolutionVSPart(outputFile_, "PhiResolutionGenVSMu", -0.002, 0.002, -0.002, 0.002);
0585   mapHisto_["PhiResolutionSimVSMu"] = new HResolutionVSPart(outputFile_, "PhiResolutionSimVSMu");
0586 
0587   // Covariances between muons kinematic quantities
0588   // ----------------------------------------------
0589   double ptMax = ptMax_;
0590 
0591   // Mass resolution
0592   // ---------------
0593   mapHisto_["MassResolution"] = new HMassResolutionVSPart(outputFile_, "MassResolution");
0594 
0595   //  mapHisto_["hResolRecoMassVSGenMassVSPt"] = new HResolutionVSPart
0596 
0597   // Mass resolution vs (pt, eta) of the muons from MC
0598   massResolutionVsPtEta_ = new HCovarianceVSxy("Mass", "Mass", 100, 0., ptMax, 60, -3, 3);
0599   // Mass resolution vs (pt, eta) of the muons from function
0600   recoPtVsgenPt_ = new TH2D("recoPtVsgenPt", "recoPtVsgenPt", 100, 0, ptMax, 100, 0, ptMax);
0601   recoPtVsgenPtEta12_ = new TH2D("recoPtVsgenPtEta12", "recoPtVsgenPtEta12", 100, 0, ptMax, 100, 0, ptMax);
0602   deltaPtOverPt_ = new TH1D("deltaPtOverPt", "deltaPtOverPt", 100, -0.1, 0.1);
0603   deltaPtOverPtForEta12_ = new TH1D("deltaPtOverPtForEta12", "deltaPtOverPtForEta12", 100, -0.1, 0.1);
0604 
0605   // Muons resolutions from resolution functions
0606   // -------------------------------------------
0607   int totBinsY = 60;
0608   mapHisto_["hFunctionResolMass"] = new HFunctionResolution(outputFile_, "hFunctionResolMass", ptMax, totBinsY);
0609   mapHisto_["hFunctionResolPt"] = new HFunctionResolution(outputFile_, "hFunctionResolPt", ptMax, totBinsY);
0610   mapHisto_["hFunctionResolCotgTheta"] =
0611       new HFunctionResolution(outputFile_, "hFunctionResolCotgTheta", ptMax, totBinsY);
0612   mapHisto_["hFunctionResolPhi"] = new HFunctionResolution(outputFile_, "hFunctionResolPhi", ptMax, totBinsY);
0613 
0614   if (readCovariances_) {
0615     // Covariances read from file. Used to compare the terms in the expression of mass error
0616     mapHisto_["ReadCovariances"] = new HCovarianceVSParts(theCovariancesRootFileName_, "Covariance");
0617 
0618     // Variances
0619     mapHisto_["MassResolutionPt"] = new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPt", ptMax);
0620     mapHisto_["MassResolutionCotgTheta"] =
0621         new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassCotgTheta", ptMax);
0622     mapHisto_["MassResolutionPhi"] = new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPhi", ptMax);
0623     // Covariances
0624     mapHisto_["MassResolutionPt-CotgTheta"] =
0625         new HFunctionResolution(outputFile_, "functionResolMassPt-CotgTheta", ptMax, totBinsY);
0626     mapHisto_["MassResolutionPt-Phi"] =
0627         new HFunctionResolution(outputFile_, "functionResolMassPt-Phi", ptMax, totBinsY);
0628     mapHisto_["MassResolutionCotgTheta-Phi"] =
0629         new HFunctionResolution(outputFile_, "functionResolMassCotgTheta-Phi", ptMax, totBinsY);
0630     mapHisto_["MassResolutionPt1-Pt2"] =
0631         new HFunctionResolution(outputFile_, "functionResolMassPt1-Pt2", ptMax, totBinsY);
0632     mapHisto_["MassResolutionCotgTheta1-CotgTheta2"] =
0633         new HFunctionResolution(outputFile_, "functionResolMassCotgTheta1-CotgTheta2", ptMax, totBinsY);
0634     mapHisto_["MassResolutionPhi1-Phi2"] =
0635         new HFunctionResolution(outputFile_, "functionResolMassPhi1-Phi2", ptMax, totBinsY);
0636     mapHisto_["MassResolutionPt12-CotgTheta21"] =
0637         new HFunctionResolution(outputFile_, "functionResolMassPt12-CotgTheta21", ptMax, totBinsY);
0638     mapHisto_["MassResolutionPt12-Phi21"] =
0639         new HFunctionResolution(outputFile_, "functionResolMassPt12-Phi21", ptMax, totBinsY);
0640     mapHisto_["MassResolutionCotgTheta12-Phi21"] =
0641         new HFunctionResolution(outputFile_, "functionResolMassCotgTheta12-Phi21", ptMax, totBinsY);
0642 
0643     mapHisto_["sigmaPtFromVariance"] = new HFunctionResolution(outputFile_, "sigmaPtFromVariance", ptMax, totBinsY);
0644     mapHisto_["sigmaCotgThetaFromVariance"] =
0645         new HFunctionResolution(outputFile_, "sigmaCotgThetaFromVariance", ptMax, totBinsY);
0646     mapHisto_["sigmaPhiFromVariance"] = new HFunctionResolution(outputFile_, "sigmaPhiFromVariance", ptMax, totBinsY);
0647 
0648     // Derivatives
0649     mapHisto_["DerivativePt"] = new HFunctionResolution(outputFile_, "derivativePt", ptMax);
0650     mapHisto_["DerivativeCotgTheta"] = new HFunctionResolution(outputFile_, "derivativeCotgTheta", ptMax);
0651     mapHisto_["DerivativePhi"] = new HFunctionResolution(outputFile_, "derivativePhi", ptMax);
0652 
0653     // Pt term from function
0654     mapHisto_["MassResolutionPtFromFunction"] =
0655         new HFunctionResolutionVarianceCheck(outputFile_, "functionResolMassPtFromFunction", ptMax);
0656 
0657     mapHisto_["FullMassResolution"] = new HFunctionResolution(outputFile_, "fullMassResolution", ptMax);
0658     mapHisto_["MassRes1"] = new HFunctionResolution(outputFile_, "massRes1", ptMax);
0659     mapHisto_["MassRes2"] = new HFunctionResolution(outputFile_, "massRes2", ptMax);
0660     mapHisto_["MassRes3"] = new HFunctionResolution(outputFile_, "massRes3", ptMax);
0661     mapHisto_["MassRes4"] = new HFunctionResolution(outputFile_, "massRes4", ptMax);
0662     mapHisto_["MassRes5"] = new HFunctionResolution(outputFile_, "massRes5", ptMax);
0663     mapHisto_["MassRes6"] = new HFunctionResolution(outputFile_, "massRes6", ptMax);
0664     mapHisto_["MassResPtAndPt12"] = new HFunctionResolution(outputFile_, "massResPtAndPt12", ptMax);
0665   } else {
0666     mapHisto_["Covariances"] = new HCovarianceVSParts(outputFile_, "Covariance", ptMax);
0667   }
0668 }
0669 
0670 void ResolutionAnalyzer::writeHistoMap() {
0671   for (std::map<std::string, Histograms*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
0672        histo++) {
0673     (*histo).second->Write();
0674   }
0675   outputFile_->cd();
0676   massResolutionVsPtEta_->Write();
0677   recoPtVsgenPt_->Write();
0678   recoPtVsgenPtEta12_->Write();
0679   deltaPtOverPt_->Write();
0680   deltaPtOverPtForEta12_->Write();
0681 }
0682 
0683 bool ResolutionAnalyzer::checkDeltaR(const reco::Particle::LorentzVector& genMu,
0684                                      const reco::Particle::LorentzVector& recMu) {
0685   //first is always mu-, second is always mu+
0686   double deltaR =
0687       sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
0688            ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
0689   if (deltaR < 0.01)
0690     return true;
0691   else if (debug_ > 0)
0692     std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
0693               << " DOES NOT MATCH with generated muon from resonance: " << std::endl
0694               << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
0695   return false;
0696 }
0697 
0698 //define this as a plug-in
0699 DEFINE_FWK_MODULE(ResolutionAnalyzer);
[ 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!