Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:32

0001 /** \class MuScleFit
0002  *  Analyzer of the Global muon tracks
0003 */
0004 
0005 //  \class MuScleFit
0006 //  Fitter of momentum scale and resolution from resonance decays to muon track pairs
0007 //
0008 //  \author R. Bellan, C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova
0009 // revised S. Casasso, E. Migliore - UniTo & INFN Torino
0010 //
0011 //  Recent additions:
0012 //  - several parameters allow a more flexible use, tests, and control handles
0013 //    for the likelihood. In particular, a set of integers controls the order
0014 //    with which parameters are released; another controls which parameters are
0015 //    fixed. A function allows to smear momenta and angles of muons from the
0016 //    resonance before any correction, using a set of random numbers generated
0017 //    once and for all in the constructor (this way the smearing remains the same
0018 //    for any given muon no matter how many times one loops and what corrections
0019 //    one applies).
0020 //    For a correct use of these flags, please see the function minimizeLikelihood() in
0021 //    MuScleFitUtils.cc
0022 //  - the fit now allows to extract resolution functions simultaneously with the
0023 //    momentum scale. So far a simple parametrization
0024 //    of muon momentum resolution and angle resolution has been implemented, but
0025 //    extensions are straightforward.
0026 //  - It is however advisable to fit separately resolution and scale. The suggested
0027 //    course of action is:
0028 //    1) fit the scale with a simple parametrization
0029 //    2) check results, fit with more complicated forms
0030 //    3) verify which is a sufficiently accurate description of the data
0031 //    4) fix scale parameters and fit for the resolution
0032 //    5) go back to fitting the scale with resolution parameters fixed to fitted values
0033 //  - Also note that resolution fits may fail to converge due to instability of the
0034 //    probability distribution to the limit of large widths. Improvements here are
0035 //    advisable.
0036 //  - The treatment of signal windows in the Y region
0037 //    has to be refined because of overlaps. More work is needed here, assigning a different
0038 //    weight to different hypothesis of the resonance producing a given mass, if there are
0039 //    multiple candidates.
0040 //  - Also, larger windows are to be allowed for fits to SA muons.
0041 //  - File Probs_1000.root contains the probability distribution of lorentzians convoluted
0042 //    with gaussian smearing functions, for the six resonances. A 1000x1000 grid
0043 //    in mass,sigma has been computed (using root macro Probs.C).
0044 //    A wider interval of masses for each resonance should be computed, to be used for standalone muons
0045 //
0046 //
0047 //  Notes on additions, TD 31/3/08
0048 //
0049 //  - background model: at least a couple of different models, with two parameters,
0050 //    should be included in the fitting procedure such that the function massprob(),
0051 //    which produces the probability in the likelihood computation, incorporates the
0052 //    probability that the event is from background. That is, if the fitting function
0053 //    knows the shape of the mass spectrum ( say, a falling exponential plus a gaussian
0054 //    signal) it becomes possible to fit the scale together with the background shape
0055 //    and normalization parameters. Of course, one should do one thing at a time: first
0056 //    a scale fit, then a shape fit with scale parameters fixed, and then a combined
0057 //    fit. Convergence issues should be handled case by case.
0058 //  - The correct implementation of the above idea requires a reorganization of pass
0059 //    parameters (in the cfg) and fit parameters. The user has to be able to smear,
0060 //    bias, fix parameters, choose scale fitting functions, resolution fitting functions,
0061 //    and background functions. It should be possible to separate the fit functions from
0062 //    the biasing ones, which would allow a more thorough testing.
0063 //  - all the above can be obtained by making the .cfg instructions heavier. Since this
0064 //    is a routine operated by experts only, it is a sensible choice.
0065 //  - One should thus envision the following:
0066 //      1) a set of parameters controlling the biasing function: parBias()
0067 //      2) a set of parameters controlling the smearing function: parSmear()
0068 //      3) a set of parameters to define resolution modeling and initial values: parResol()
0069 //      3b) parResol() gets fix and order bits by parResolFix() and parResolOrder()
0070 //      4) a set of parameters to define scale modeling and initial values: parScale()
0071 //      4b) parScale() gets fix and order bits by parScaleFix() and parScaleOrder()
0072 //      5) a set of parameters controlling the background shape and normalization: parNorm()
0073 //      5b) parNorm() gets fix and order bits by parNormFix() and parNormOrder()
0074 //    The likelihood parameters then become a vector which is dynamically composed of
0075 //    sets 3), 4), and 5): parval() = parResol()+parScale()+parNorm()
0076 //  - In order to study better the likelihood behavior it would be advisable to introduce
0077 //    some histogram filling on the last iteration of the likelihood. It is not clear
0078 //    how best to achieve that: probably the simplest way is to make a histogram filling
0079 //    function run just after the likelihood computation, such that the best value of the
0080 //    fit parameters is used.
0081 //  - The muon pair which we call our resonance must be chosen in a way which does not
0082 //    bias our likelihood: we cannot just choose the pair closest to a resonance.
0083 //
0084 //
0085 //    Notes on additions, T.Dorigo 22/12/2008
0086 //    ---------------------------------------
0087 //
0088 //  - File Probs_new_1000_CTEQ.root now contains a set of 24 additional two-dim histograms,
0089 //    defining the probability distribution of Z boson decays as a function of measured mass
0090 //    and expected sigma in 24 different bins of Z rapidity, extracted from CTEQ 6 PDF (at
0091 //    Leading Order) from the convolution in the factorization integral. See programs CTEQ.cpp
0092 //    and Fits.C.
0093 //  - The probability for Z boson events now thus depends on the measured rapidity of the dimuon
0094 //    system. All functions in file MuScleFitUtils.cc have been suitably changed.
0095 //
0096 // ----------------------------------------------------------------------------------
0097 //    Modifications by M. De Mattia 13/3/2009
0098 //    ---------------------------------------
0099 //  - The histograms map was moved to a base class (MuScleFitBase) from which this one inherits.
0100 //
0101 //    Modifications by M. De Mattia 20/7/2009
0102 //    ---------------------------------------
0103 //  - Reworked background fit based on ranges. See comments in the code for more details.
0104 // ---------------------------------------------------------------------------------------------
0105 // Base Class Headers
0106 // ------------------
0107 #include "FWCore/Framework/interface/EDLooper.h"
0108 #include "FWCore/Framework/interface/LooperFactory.h"
0109 #include "FWCore/Utilities/interface/InputTag.h"
0110 #include "FWCore/Framework/interface/Frameworkfwd.h"
0111 #include "FWCore/Framework/interface/Event.h"
0112 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0113 
0114 #include <CLHEP/Vector/LorentzVector.h>
0115 #include <memory>
0116 
0117 #include <vector>
0118 
0119 #include "FWCore/Framework/interface/EventSetup.h"
0120 #include "FWCore/Framework/interface/ESHandle.h"
0121 #include "FWCore/ServiceRegistry/interface/Service.h"
0122 
0123 #include "MuScleFitBase.h"
0124 #include "Histograms.h"
0125 #include "MuScleFitPlotter.h"
0126 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
0127 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
0128 #include "MuonAnalysis/MomentumScaleCalibration/interface/Muon.h"
0129 #include "MuonAnalysis/MomentumScaleCalibration/interface/Event.h"
0130 #include "MuScleFitMuonSelector.h"
0131 
0132 #include "DataFormats/TrackReco/interface/Track.h"
0133 #include "DataFormats/MuonReco/interface/Muon.h"
0134 #include "DataFormats/MuonReco/interface/CaloMuon.h"
0135 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0136 
0137 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0138 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0139 
0140 #include "DataFormats/PatCandidates/interface/Muon.h"
0141 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0142 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0143 #include "DataFormats/Common/interface/TriggerResults.h"
0144 #include "FWCore/Common/interface/TriggerNames.h"
0145 
0146 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0147 
0148 #include "HepPDT/defs.h"
0149 #include "HepPDT/TableBuilder.hh"
0150 #include "HepPDT/ParticleDataTable.hh"
0151 
0152 #include "HepMC/GenParticle.h"
0153 #include "HepMC/GenEvent.h"
0154 
0155 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0156 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0157 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0158 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0159 
0160 #include "DataFormats/VertexReco/interface/Vertex.h"
0161 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0162 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0163 
0164 #include "TFile.h"
0165 #include "TTree.h"
0166 #include "TMinuit.h"
0167 
0168 // To use callgrind for code profiling uncomment also the following define.
0169 // #define USE_CALLGRIND
0170 
0171 #ifdef USE_CALLGRIND
0172 #include "valgrind/callgrind.h"
0173 #endif
0174 
0175 // To read likelihood distributions from the database.
0176 //#include "CondFormats/RecoMuonObjects/interface/MuScleFitLikelihoodPdf.h"
0177 //#include "CondFormats/DataRecord/interface/MuScleFitLikelihoodPdfRcd.h"
0178 
0179 namespace edm {
0180   class ParameterSet;
0181   class Event;
0182   class EventSetup;
0183 }  // namespace edm
0184 
0185 class MuScleFit : public edm::EDLooper, MuScleFitBase {
0186 public:
0187   // Constructor
0188   // -----------
0189   MuScleFit(const edm::ParameterSet& pset);
0190 
0191   // Destructor
0192   // ----------
0193   ~MuScleFit() override;
0194 
0195   // Operations
0196   // ----------
0197   void beginOfJobInConstructor();
0198   // void beginOfJob( const edm::EventSetup& eventSetup );
0199   // virtual void beginOfJob();
0200   void endOfJob() override;
0201 
0202   void startingNewLoop(unsigned int iLoop) override;
0203 
0204   edm::EDLooper::Status endOfLoop(const edm::EventSetup& eventSetup, unsigned int iLoop) override;
0205   virtual void endOfFastLoop(const unsigned int iLoop);
0206 
0207   edm::EDLooper::Status duringLoop(const edm::Event& event, const edm::EventSetup& eventSetup) override;
0208   /**
0209    * This method performs all needed operations on the muon pair. It reads the muons from SavedPair and uses the iev
0210    * counter to keep track of the event number. The iev is incremented internally and reset to 0 in startingNewLoop.
0211    */
0212   virtual void duringFastLoop();
0213 
0214   template <typename T>
0215   std::vector<MuScleFitMuon> fillMuonCollection(const std::vector<T>& tracks);
0216 
0217 private:
0218 protected:
0219   /**
0220    * Selects the muon pairs and fills the SavedPair and (if needed) the genPair vector.
0221    * This version reads the events from the edm root file and performs a selection of the muons according to the parameters in the cfg.
0222    */
0223   void selectMuons(const edm::Event& event);
0224   /**
0225    * Selects the muon pairs and fills the SavedPair and (if needed) the genPair vector.
0226    * This version reads the events from a tree in the file specified in the cfg. The tree only contains one muon pair per event. This
0227    * means that no selection is performed and we use preselected muons.
0228    */
0229   void selectMuons(const int maxEvents, const TString& treeFileName);
0230 
0231   /// Template method used to fill the track collection starting from reco::muons or pat::muons
0232   template <typename T>
0233   void takeSelectedMuonType(const T& muon, std::vector<reco::Track>& tracks);
0234   /// Function for onia selections
0235   bool selGlobalMuon(const pat::Muon* aMuon);
0236   bool selTrackerMuon(const pat::Muon* aMuon);
0237 
0238   /// Check if two lorentzVector are near in deltaR
0239   bool checkDeltaR(reco::Particle::LorentzVector& genMu, reco::Particle::LorentzVector& recMu);
0240   /// Fill the reco vs gen and reco vs sim comparison histograms
0241   void fillComparisonHistograms(const reco::Particle::LorentzVector& genMu,
0242                                 const reco::Particle::LorentzVector& recoMu,
0243                                 const std::string& inputName,
0244                                 const int charge);
0245 
0246   /// Apply the smearing if needed using the function in MuScleFitUtils
0247   void applySmearing(reco::Particle::LorentzVector& mu);
0248   /// Apply the bias if needed using the function in MuScleFitUtils
0249   void applyBias(reco::Particle::LorentzVector& mu, const int charge);
0250 
0251   /**
0252    * Simple method to check parameters consistency. It aborts the job if the parameters
0253    * are not consistent.
0254    */
0255   void checkParameters();
0256 
0257   MuonServiceProxy* theService;
0258 
0259   // Counters
0260   // --------
0261   int numberOfSimTracks;
0262   int numberOfSimMuons;
0263   int numberOfSimVertices;
0264   int numberOfEwkZ;
0265 
0266   bool ifHepMC;
0267   bool ifGenPart;
0268 
0269   // Constants
0270   // ---------
0271   double minResMass_hwindow[6];
0272   double maxResMass_hwindow[6];
0273 
0274   // Total number of loops
0275   // ---------------------
0276   unsigned int maxLoopNumber;
0277   unsigned int loopCounter;
0278 
0279   bool fastLoop;
0280 
0281   MuScleFitPlotter* plotter;
0282 
0283   // The reconstructed muon 4-momenta to be put in the tree
0284   // ------------------------------------------------------
0285   reco::Particle::LorentzVector recMu1, recMu2;
0286   MuScleFitMuon recMuScleMu1, recMuScleMu2;
0287   int iev;
0288   int totalEvents_;
0289 
0290   bool compareToSimTracks_;
0291   edm::InputTag simTracksCollection_;
0292   bool PATmuons_;
0293   std::string genParticlesName_;
0294 
0295   // Input Root Tree file name. If empty events are read from the edm root file.
0296   std::string inputRootTreeFileName_;
0297   // Output Root Tree file name. If not empty events are dumped to this file at the end of the last iteration.
0298   std::string outputRootTreeFileName_;
0299   // Maximum number of events from root tree. It works in the same way as the maxEvents to configure a input source.
0300   int maxEventsFromRootTree_;
0301 
0302   std::string triggerResultsLabel_;
0303   std::string triggerResultsProcess_;
0304   std::vector<std::string> triggerPath_;
0305   bool negateTrigger_;
0306   bool saveAllToTree_;
0307 
0308   // input collections for PU related infos
0309   edm::InputTag puInfoSrc_;
0310   edm::InputTag vertexSrc_;
0311 
0312   std::unique_ptr<MuScleFitMuonSelector> muonSelector_;
0313 };
0314 
0315 template <typename T>
0316 std::vector<MuScleFitMuon> MuScleFit::fillMuonCollection(const std::vector<T>& tracks) {
0317   std::vector<MuScleFitMuon> muons;
0318   typename std::vector<T>::const_iterator track;
0319   for (track = tracks.begin(); track != tracks.end(); ++track) {
0320     reco::Particle::LorentzVector mu;
0321     mu = reco::Particle::LorentzVector(
0322         track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
0323     // Apply smearing if needed, and then bias
0324     // ---------------------------------------
0325     MuScleFitUtils::goodmuon++;
0326     if (debug_ > 0)
0327       std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon << ": initial value   Pt = " << mu.Pt()
0328                 << std::endl;
0329 
0330     applySmearing(mu);
0331     applyBias(mu, track->charge());
0332     if (debug_ > 0)
0333       std::cout << "track charge: " << track->charge() << std::endl;
0334 
0335     Double_t hitsTk = track->innerTrack()->hitPattern().numberOfValidTrackerHits();
0336     Double_t hitsMuon = track->innerTrack()->hitPattern().numberOfValidMuonHits();
0337     Double_t ptError = track->innerTrack()->ptError();
0338     MuScleFitMuon muon(mu, track->charge(), ptError, hitsTk, hitsMuon, false);
0339     if (debug_ > 0) {
0340       std::cout << "[MuScleFit::fillMuonCollection]" << std::endl;
0341       std::cout << "  muon = " << muon << std::endl;
0342     }
0343 
0344     // Store modified muon
0345     // -------------------
0346     muons.push_back(muon);
0347   }
0348   return muons;
0349 }
0350 
0351 template <typename T>
0352 void MuScleFit::takeSelectedMuonType(const T& muon, std::vector<reco::Track>& tracks) {
0353   // std::cout<<"muon "<<muon->isGlobalMuon()<<muon->isStandAloneMuon()<<muon->isTrackerMuon()<<std::endl;
0354   //NNBB: one muon can be of many kinds at once but with the theMuonType_ we are sure
0355   // to avoid double counting of the same muon
0356   if (muon->isGlobalMuon() && theMuonType_ == 1)
0357     tracks.push_back(*(muon->globalTrack()));
0358   else if (muon->isStandAloneMuon() && theMuonType_ == 2)
0359     tracks.push_back(*(muon->outerTrack()));
0360   else if (muon->isTrackerMuon() && theMuonType_ == 3)
0361     tracks.push_back(*(muon->innerTrack()));
0362 
0363   else if (theMuonType_ == 10 && !(muon->isStandAloneMuon()))  //particular case!!
0364     tracks.push_back(*(muon->innerTrack()));
0365   else if (theMuonType_ == 11 && muon->isGlobalMuon())
0366     tracks.push_back(*(muon->innerTrack()));
0367   else if (theMuonType_ == 13 && muon->isTrackerMuon())
0368     tracks.push_back(*(muon->innerTrack()));
0369 }
0370 
0371 // Constructor
0372 // -----------
0373 MuScleFit::MuScleFit(const edm::ParameterSet& pset) : MuScleFitBase(pset), totalEvents_(0) {
0374   MuScleFitUtils::debug = debug_;
0375   if (debug_ > 0)
0376     std::cout << "[MuScleFit]: Constructor" << std::endl;
0377 
0378   if ((theMuonType_ < -4 || theMuonType_ > 5) && theMuonType_ < 10) {
0379     std::cout << "[MuScleFit]: Unknown muon type! Aborting." << std::endl;
0380     abort();
0381   }
0382 
0383   loopCounter = 0;
0384 
0385   // Boundaries for h-function computation (to be improved!)
0386   // -------------------------------------------------------
0387   minResMass_hwindow[0] = 71.1876;  // 76.;
0388   maxResMass_hwindow[0] = 111.188;  // 106.;
0389   minResMass_hwindow[1] = 10.15;
0390   maxResMass_hwindow[1] = 10.55;
0391   minResMass_hwindow[2] = 9.8;
0392   maxResMass_hwindow[2] = 10.2;
0393   minResMass_hwindow[3] = 9.25;
0394   maxResMass_hwindow[3] = 9.65;
0395   minResMass_hwindow[4] = 3.58;
0396   maxResMass_hwindow[4] = 3.78;
0397   minResMass_hwindow[5] = 3.0;
0398   maxResMass_hwindow[5] = 3.2;
0399 
0400   // Max number of loops (if > 2 then try to minimize likelihood more than once)
0401   // ---------------------------------------------------------------------------
0402   maxLoopNumber = pset.getUntrackedParameter<int>("maxLoopNumber", 2);
0403   fastLoop = pset.getUntrackedParameter<bool>("FastLoop", true);
0404 
0405   // Selection of fits according to loop
0406   MuScleFitUtils::doResolFit = pset.getParameter<std::vector<int> >("doResolFit");
0407   MuScleFitUtils::doScaleFit = pset.getParameter<std::vector<int> >("doScaleFit");
0408   MuScleFitUtils::doCrossSectionFit = pset.getParameter<std::vector<int> >("doCrossSectionFit");
0409   MuScleFitUtils::doBackgroundFit = pset.getParameter<std::vector<int> >("doBackgroundFit");
0410 
0411   // Bias and smear types
0412   // --------------------
0413   int biasType = pset.getParameter<int>("BiasType");
0414   MuScleFitUtils::BiasType = biasType;
0415   // No error, the scale functions are used also for the bias
0416   MuScleFitUtils::biasFunction = scaleFunctionVecService(biasType);
0417   int smearType = pset.getParameter<int>("SmearType");
0418   MuScleFitUtils::SmearType = smearType;
0419   MuScleFitUtils::smearFunction = smearFunctionService(smearType);
0420 
0421   // Fit types
0422   // ---------
0423   int resolFitType = pset.getParameter<int>("ResolFitType");
0424   MuScleFitUtils::ResolFitType = resolFitType;
0425   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType);
0426   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolFitType);
0427   int scaleType = pset.getParameter<int>("ScaleFitType");
0428   MuScleFitUtils::ScaleFitType = scaleType;
0429   MuScleFitUtils::scaleFunction = scaleFunctionService(scaleType);
0430   MuScleFitUtils::scaleFunctionForVec = scaleFunctionVecService(scaleType);
0431 
0432   // Initial parameters values
0433   // -------------------------
0434   MuScleFitUtils::parBias = pset.getParameter<std::vector<double> >("parBias");
0435   MuScleFitUtils::parSmear = pset.getParameter<std::vector<double> >("parSmear");
0436   MuScleFitUtils::parResol = pset.getParameter<std::vector<double> >("parResol");
0437   MuScleFitUtils::parResolStep =
0438       pset.getUntrackedParameter<std::vector<double> >("parResolStep", std::vector<double>());
0439   MuScleFitUtils::parResolMin = pset.getUntrackedParameter<std::vector<double> >("parResolMin", std::vector<double>());
0440   MuScleFitUtils::parResolMax = pset.getUntrackedParameter<std::vector<double> >("parResolMax", std::vector<double>());
0441   MuScleFitUtils::parScale = pset.getParameter<std::vector<double> >("parScale");
0442   MuScleFitUtils::parScaleStep =
0443       pset.getUntrackedParameter<std::vector<double> >("parScaleStep", std::vector<double>());
0444   MuScleFitUtils::parScaleMin = pset.getUntrackedParameter<std::vector<double> >("parScaleMin", std::vector<double>());
0445   MuScleFitUtils::parScaleMax = pset.getUntrackedParameter<std::vector<double> >("parScaleMax", std::vector<double>());
0446   MuScleFitUtils::parCrossSection = pset.getParameter<std::vector<double> >("parCrossSection");
0447   MuScleFitUtils::parBgr = pset.getParameter<std::vector<double> >("parBgr");
0448   MuScleFitUtils::parResolFix = pset.getParameter<std::vector<int> >("parResolFix");
0449   MuScleFitUtils::parScaleFix = pset.getParameter<std::vector<int> >("parScaleFix");
0450   MuScleFitUtils::parCrossSectionFix = pset.getParameter<std::vector<int> >("parCrossSectionFix");
0451   MuScleFitUtils::parBgrFix = pset.getParameter<std::vector<int> >("parBgrFix");
0452   MuScleFitUtils::parResolOrder = pset.getParameter<std::vector<int> >("parResolOrder");
0453   MuScleFitUtils::parScaleOrder = pset.getParameter<std::vector<int> >("parScaleOrder");
0454   MuScleFitUtils::parCrossSectionOrder = pset.getParameter<std::vector<int> >("parCrossSectionOrder");
0455   MuScleFitUtils::parBgrOrder = pset.getParameter<std::vector<int> >("parBgrOrder");
0456 
0457   MuScleFitUtils::resfind = pset.getParameter<std::vector<int> >("resfind");
0458   MuScleFitUtils::FitStrategy = pset.getParameter<int>("FitStrategy");
0459 
0460   // Option to skip unnecessary stuff
0461   // --------------------------------
0462   MuScleFitUtils::speedup = pset.getParameter<bool>("speedup");
0463 
0464   // Option to skip simTracks comparison
0465   compareToSimTracks_ = pset.getParameter<bool>("compareToSimTracks");
0466   simTracksCollection_ = pset.getUntrackedParameter<edm::InputTag>("SimTracksCollection", edm::InputTag("g4SimHits"));
0467 
0468   triggerResultsLabel_ = pset.getUntrackedParameter<std::string>("TriggerResultsLabel");
0469   triggerResultsProcess_ = pset.getUntrackedParameter<std::string>("TriggerResultsProcess");
0470   triggerPath_ = pset.getUntrackedParameter<std::vector<std::string> >("TriggerPath");
0471   negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
0472   saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
0473 
0474   // input collections for PU related infos
0475   puInfoSrc_ = pset.getUntrackedParameter<edm::InputTag>("PileUpSummaryInfo");
0476   vertexSrc_ = pset.getUntrackedParameter<edm::InputTag>("PrimaryVertexCollection");
0477 
0478   PATmuons_ = pset.getUntrackedParameter<bool>("PATmuons", false);
0479   genParticlesName_ = pset.getUntrackedParameter<std::string>("GenParticlesName", "genParticles");
0480 
0481   // Use the probability file or not. If not it will perform a simpler selection taking the muon pair with
0482   // invariant mass closer to the pdf value and will crash if some fit is attempted.
0483   MuScleFitUtils::useProbsFile_ = pset.getUntrackedParameter<bool>("UseProbsFile", true);
0484 
0485   // This must be set to true if using events generated with Sherpa
0486   MuScleFitUtils::sherpa_ = pset.getUntrackedParameter<bool>("Sherpa", false);
0487 
0488   MuScleFitUtils::rapidityBinsForZ_ = pset.getUntrackedParameter<bool>("RapidityBinsForZ", true);
0489 
0490   // Set the cuts on muons to be used in the fit
0491   MuScleFitUtils::separateRanges_ = pset.getUntrackedParameter<bool>("SeparateRanges", true);
0492   MuScleFitUtils::maxMuonPt_ = pset.getUntrackedParameter<double>("MaxMuonPt", 100000000.);
0493   MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
0494   MuScleFitUtils::minMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MinMuonEtaFirstRange", -6.);
0495   MuScleFitUtils::maxMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaFirstRange", 6.);
0496   MuScleFitUtils::minMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MinMuonEtaSecondRange", -100.);
0497   MuScleFitUtils::maxMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaSecondRange", 100.);
0498   MuScleFitUtils::deltaPhiMinCut_ = pset.getUntrackedParameter<double>("DeltaPhiMinCut", -100.);
0499   MuScleFitUtils::deltaPhiMaxCut_ = pset.getUntrackedParameter<double>("DeltaPhiMaxCut", 100.);
0500 
0501   MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("DebugMassResol", false);
0502   // MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents;
0503 
0504   // Check for parameters consistency
0505   // it will abort in case of errors.
0506   checkParameters();
0507 
0508   // Generate array of gaussian-distributed numbers for smearing
0509   // -----------------------------------------------------------
0510   if (MuScleFitUtils::SmearType > 0) {
0511     std::cout << "[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
0512     TF1 G("G", "[0]*exp(-0.5*pow(x,2))", -5., 5.);
0513     double norm = 1 / sqrt(2 * TMath::Pi());
0514     G.SetParameter(0, norm);
0515     for (int i = 0; i < 10000; i++) {
0516       for (int j = 0; j < 7; j++) {
0517         MuScleFitUtils::x[j][i] = G.GetRandom();
0518       }
0519     }
0520   }
0521   MuScleFitUtils::goodmuon = 0;
0522 
0523   if (theMuonType_ > 0 && theMuonType_ < 4) {
0524     MuScleFitUtils::MuonTypeForCheckMassWindow = theMuonType_ - 1;
0525     MuScleFitUtils::MuonType = theMuonType_ - 1;
0526   } else if (theMuonType_ == 0 || theMuonType_ == 4 || theMuonType_ == 5 || theMuonType_ >= 10 || theMuonType_ == -1 ||
0527              theMuonType_ == -2 || theMuonType_ == -3 || theMuonType_ == -4) {
0528     MuScleFitUtils::MuonTypeForCheckMassWindow = 2;
0529     MuScleFitUtils::MuonType = 2;
0530   } else {
0531     std::cout << "Wrong muon type " << theMuonType_ << std::endl;
0532     exit(1);
0533   }
0534 
0535   // When using standalone muons switch to the single Z pdf
0536   if (theMuonType_ == 2) {
0537     MuScleFitUtils::rapidityBinsForZ_ = false;
0538   }
0539 
0540   // Initialize ResMaxSigma And ResHalfWidth - 0 = global, 1 = SM, 2 = tracker
0541   // -------------------------------------------------------------------------
0542   MuScleFitUtils::massWindowHalfWidth[0][0] = 20.;
0543   MuScleFitUtils::massWindowHalfWidth[0][1] = 0.35;
0544   MuScleFitUtils::massWindowHalfWidth[0][2] = 0.35;
0545   MuScleFitUtils::massWindowHalfWidth[0][3] = 0.35;
0546   MuScleFitUtils::massWindowHalfWidth[0][4] = 0.2;
0547   MuScleFitUtils::massWindowHalfWidth[0][5] = 0.2;
0548   MuScleFitUtils::massWindowHalfWidth[1][0] = 20.;
0549   MuScleFitUtils::massWindowHalfWidth[1][1] = 0.35;
0550   MuScleFitUtils::massWindowHalfWidth[1][2] = 0.35;
0551   MuScleFitUtils::massWindowHalfWidth[1][3] = 0.35;
0552   MuScleFitUtils::massWindowHalfWidth[1][4] = 0.2;
0553   MuScleFitUtils::massWindowHalfWidth[1][5] = 0.2;
0554   MuScleFitUtils::massWindowHalfWidth[2][0] = 20.;
0555   MuScleFitUtils::massWindowHalfWidth[2][1] = 0.35;
0556   MuScleFitUtils::massWindowHalfWidth[2][2] = 0.35;
0557   MuScleFitUtils::massWindowHalfWidth[2][3] = 0.35;
0558   MuScleFitUtils::massWindowHalfWidth[2][4] = 0.2;
0559   MuScleFitUtils::massWindowHalfWidth[2][5] = 0.2;
0560 
0561   muonSelector_ = std::make_unique<MuScleFitMuonSelector>(theMuonLabel_,
0562                                                           theMuonType_,
0563                                                           PATmuons_,
0564                                                           MuScleFitUtils::resfind,
0565                                                           MuScleFitUtils::speedup,
0566                                                           genParticlesName_,
0567                                                           compareToSimTracks_,
0568                                                           simTracksCollection_,
0569                                                           MuScleFitUtils::sherpa_,
0570                                                           debug_);
0571 
0572   MuScleFitUtils::backgroundHandler =
0573       new BackgroundHandler(pset.getParameter<std::vector<int> >("BgrFitType"),
0574                             pset.getParameter<std::vector<double> >("LeftWindowBorder"),
0575                             pset.getParameter<std::vector<double> >("RightWindowBorder"),
0576                             MuScleFitUtils::ResMass,
0577                             MuScleFitUtils::massWindowHalfWidth[MuScleFitUtils::MuonTypeForCheckMassWindow]);
0578 
0579   MuScleFitUtils::crossSectionHandler =
0580       new CrossSectionHandler(MuScleFitUtils::parCrossSection, MuScleFitUtils::resfind);
0581 
0582   // Build cross section scale factors
0583   // MuScleFitUtils::resfind
0584 
0585   MuScleFitUtils::normalizeLikelihoodByEventNumber_ =
0586       pset.getUntrackedParameter<bool>("NormalizeLikelihoodByEventNumber", true);
0587   if (debug_ > 0)
0588     std::cout << "End of MuScleFit constructor" << std::endl;
0589 
0590   inputRootTreeFileName_ = pset.getParameter<std::string>("InputRootTreeFileName");
0591   outputRootTreeFileName_ = pset.getParameter<std::string>("OutputRootTreeFileName");
0592   maxEventsFromRootTree_ = pset.getParameter<int>("MaxEventsFromRootTree");
0593 
0594   MuScleFitUtils::startWithSimplex_ = pset.getParameter<bool>("StartWithSimplex");
0595   MuScleFitUtils::computeMinosErrors_ = pset.getParameter<bool>("ComputeMinosErrors");
0596   MuScleFitUtils::minimumShapePlots_ = pset.getParameter<bool>("MinimumShapePlots");
0597 
0598   beginOfJobInConstructor();
0599 }
0600 
0601 // Destructor
0602 // ----------
0603 MuScleFit::~MuScleFit() {
0604   if (debug_ > 0)
0605     std::cout << "[MuScleFit]: Destructor" << std::endl;
0606   std::cout << "Total number of analyzed events = " << totalEvents_ << std::endl;
0607 
0608   if (!(outputRootTreeFileName_.empty())) {
0609     // Save the events to a root tree unless we are reading from the edm root file and the SavedPair size is different from the totalEvents_
0610     if (!(inputRootTreeFileName_.empty() && (int(MuScleFitUtils::SavedPair.size()) != totalEvents_))) {
0611       std::cout << "Saving muon pairs to root tree" << std::endl;
0612       RootTreeHandler rootTreeHandler;
0613       if (MuScleFitUtils::speedup) {
0614         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, 0, saveAllToTree_);
0615         if (debug_ > 0) {
0616           std::vector<MuonPair>::const_iterator it = muonPairs_.begin();
0617           std::cout << "[MuScleFit::~MuScleFit] (Destructor)" << std::endl;
0618           for (; it < muonPairs_.end(); ++it) {
0619             std::cout << "  Debugging pairs that are going to be written to file" << std::endl;
0620             std::cout << "  muon1 = " << it->mu1 << std::endl;
0621             std::cout << "  muon2 = " << it->mu2 << std::endl;
0622           }
0623         }
0624         rootTreeHandler.writeTree(outputRootTreeFileName_, &(muonPairs_), theMuonType_, nullptr, saveAllToTree_);
0625       } else {
0626         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, &(MuScleFitUtils::genPair), saveAllToTree_ );
0627         rootTreeHandler.writeTree(
0628             outputRootTreeFileName_, &(muonPairs_), theMuonType_, &(genMuonPairs_), saveAllToTree_);
0629       }
0630     } else {
0631       std::cout << "ERROR: events in the vector = " << MuScleFitUtils::SavedPair.size()
0632                 << " != totalEvents = " << totalEvents_ << std::endl;
0633     }
0634   }
0635 }
0636 
0637 // Begin job
0638 // ---------
0639 void MuScleFit::beginOfJobInConstructor()
0640 // void MuScleFit::beginOfJob ()
0641 // void MuScleFit::beginOfJob( const edm::EventSetup& eventSetup )
0642 {
0643   if (debug_ > 0)
0644     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0645   //if(maxLoopNumber>1)
0646   if (MuScleFitUtils::useProbsFile_) {
0647     readProbabilityDistributionsFromFile();
0648   }
0649 
0650   if (debug_ > 0)
0651     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0652 
0653   // Create the root file
0654   // --------------------
0655   for (unsigned int i = 0; i < (maxLoopNumber); i++) {
0656     std::stringstream ss;
0657     ss << i;
0658     std::string rootFileName = ss.str() + "_" + theRootFileName_;
0659     if (theCompressionSettings_ > -1) {
0660       theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE", "", theCompressionSettings_));
0661     } else {
0662       theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE"));
0663     }
0664   }
0665   if (debug_ > 0)
0666     std::cout << "[MuScleFit]: Root file created" << std::endl;
0667 
0668   std::cout << "creating plotter" << std::endl;
0669   plotter = new MuScleFitPlotter(theGenInfoRootFileName_);
0670   plotter->debug = debug_;
0671 }
0672 
0673 // End of job method
0674 // -----------------
0675 void MuScleFit::endOfJob() {
0676   if (debug_ > 0)
0677     std::cout << "[MuScleFit]: endOfJob" << std::endl;
0678 }
0679 
0680 // New loop
0681 // --------
0682 void MuScleFit::startingNewLoop(unsigned int iLoop) {
0683   if (debug_ > 0)
0684     std::cout << "[MuScleFit]: Starting loop # " << iLoop << std::endl;
0685 
0686   // Number of muons used
0687   // --------------------
0688   MuScleFitUtils::goodmuon = 0;
0689 
0690   // Counters for problem std::cout-ing
0691   // -----------------------------
0692   MuScleFitUtils::counter_resprob = 0;
0693 
0694   // Create the root file
0695   // --------------------
0696   fillHistoMap(theFiles_[iLoop], iLoop);
0697 
0698   loopCounter = iLoop;
0699   MuScleFitUtils::loopCounter = loopCounter;
0700 
0701   iev = 0;
0702   MuScleFitUtils::iev_ = 0;
0703 
0704   MuScleFitUtils::oldNormalization_ = 0;
0705 }
0706 
0707 // End of loop routine
0708 // -------------------
0709 edm::EDLooper::Status MuScleFit::endOfLoop(const edm::EventSetup& eventSetup, unsigned int iLoop) {
0710   unsigned int iFastLoop = 1;
0711 
0712   // Read the events from the root tree if requested
0713   if (!(inputRootTreeFileName_.empty())) {
0714     selectMuons(maxEventsFromRootTree_, inputRootTreeFileName_);
0715     // When reading from local file all the loops are done here
0716     totalEvents_ = MuScleFitUtils::SavedPair.size();
0717     iFastLoop = 0;
0718   } else {
0719     endOfFastLoop(iLoop);
0720   }
0721 
0722   // If a fastLoop is required we do all the remaining iterations here
0723   if (fastLoop == true) {
0724     for (; iFastLoop < maxLoopNumber; ++iFastLoop) {
0725       std::cout << "Starting fast loop number " << iFastLoop << std::endl;
0726 
0727       // In the first loop is called by the framework
0728       // if( iFastLoop > 0 ) {
0729       startingNewLoop(iFastLoop);
0730       // }
0731 
0732       // std::vector<std::pair<lorentzVector,lorentzVector> >::const_iterator it = MuScleFitUtils::SavedPair.begin();
0733       // for( ; it != SavedPair.end(); ++it ) {
0734       while (iev < totalEvents_) {
0735         if (iev % 50000 == 0) {
0736           std::cout << "Fast looping on event number " << iev << std::endl;
0737         }
0738         // This reads muons from SavedPair using iev to keep track of the event
0739         duringFastLoop();
0740       }
0741       std::cout << "End of fast loop number " << iFastLoop << ". Ran on " << iev << " events" << std::endl;
0742       endOfFastLoop(iFastLoop);
0743     }
0744   }
0745 
0746   if (iFastLoop >= maxLoopNumber - 1) {
0747     return kStop;
0748   } else {
0749     return kContinue;
0750   }
0751 }
0752 
0753 void MuScleFit::endOfFastLoop(const unsigned int iLoop) {
0754   // std::cout<< "Inside endOfFastLoop, iLoop = " << iLoop << " and loopCounter = " << loopCounter << std::endl;
0755 
0756   if (loopCounter == 0) {
0757     // plotter->writeHistoMap();
0758     // The destructor will call the writeHistoMap after the cd to the output file
0759     delete plotter;
0760   }
0761 
0762   std::cout << "Ending loop # " << iLoop << std::endl;
0763 
0764   // Write the histos to file
0765   // ------------------------
0766   // theFiles_[iLoop]->cd();
0767   writeHistoMap(iLoop);
0768 
0769   // Likelihood minimization to compute corrections
0770   // ----------------------------------------------
0771   // theFiles_[iLoop]->cd();
0772   TDirectory* likelihoodDir = theFiles_[iLoop]->mkdir("likelihood");
0773   likelihoodDir->cd();
0774   MuScleFitUtils::minimizeLikelihood();
0775 
0776   // ATTENTION, this was put BEFORE the minimizeLikelihood. Check for problems.
0777   theFiles_[iLoop]->Close();
0778   // ATTENTION: Check that this delete does not give any problem
0779   delete theFiles_[iLoop];
0780 
0781   // Clear the histos
0782   // ----------------
0783   clearHistoMap();
0784 }
0785 
0786 // Stuff to do during loop
0787 // -----------------------
0788 edm::EDLooper::Status MuScleFit::duringLoop(const edm::Event& event, const edm::EventSetup& eventSetup) {
0789   edm::Handle<edm::TriggerResults> triggerResults;
0790   event.getByLabel(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()), triggerResults);
0791   //event.getByLabel(InputTag(triggerResultsLabel_),triggerResults);
0792   bool isFired = false;
0793 
0794   if (triggerPath_[0].empty())
0795     isFired = true;
0796   else if (triggerPath_[0] == "All") {
0797     isFired = triggerResults->accept();
0798     if (debug_ > 0)
0799       std::cout << "Trigger " << isFired << std::endl;
0800   } else {
0801     bool changed;
0802     HLTConfigProvider hltConfig;
0803     hltConfig.init(event.getRun(), eventSetup, triggerResultsProcess_, changed);
0804 
0805     const edm::TriggerNames& triggerNames = event.triggerNames(*triggerResults);
0806 
0807     for (unsigned i = 0; i < triggerNames.size(); i++) {
0808       const std::string& hltName = triggerNames.triggerName(i);
0809 
0810       // match the path in the pset with the true name of the trigger
0811       for (unsigned int ipath = 0; ipath < triggerPath_.size(); ipath++) {
0812         if (hltName.find(triggerPath_[ipath]) != std::string::npos) {
0813           unsigned int triggerIndex(hltConfig.triggerIndex(hltName));
0814 
0815           // triggerIndex must be less than the size of HLTR or you get a CMSException: _M_range_check
0816           if (triggerIndex < triggerResults->size()) {
0817             isFired = triggerResults->accept(triggerIndex);
0818             if (debug_ > 0)
0819               std::cout << triggerPath_[ipath] << " " << hltName << " " << isFired << std::endl;
0820           }
0821         }  // end if (matching the path in the pset with the true trigger name
0822       }
0823     }
0824   }
0825 
0826   if (negateTrigger_ && isFired)
0827     return kContinue;
0828   else if (!(negateTrigger_) && !isFired)
0829     return kContinue;
0830 
0831 #ifdef USE_CALLGRIND
0832   CALLGRIND_START_INSTRUMENTATION;
0833 #endif
0834 
0835   if (debug_ > 0) {
0836     std::cout << "[MuScleFit-duringLoop]: loopCounter = " << loopCounter << " Run: " << event.id().run()
0837               << " Event: " << event.id().event() << std::endl;
0838   }
0839 
0840   // On the first iteration we read the bank, otherwise we fetch the information from the muon tree
0841   // ------------------------------------ Important Note --------------------------------------- //
0842   // The fillMuonCollection method applies any smearing or bias to the muons, so we NEVER use
0843   // unbiased muons.
0844   // ----------------------------------------------------------------------------------------------
0845   if (loopCounter == 0) {
0846     if (!fastLoop || inputRootTreeFileName_.empty()) {
0847       if (debug_ > 0)
0848         std::cout << "Reading from edm event" << std::endl;
0849       selectMuons(event);
0850       duringFastLoop();
0851       ++totalEvents_;
0852     }
0853   }
0854 
0855   return kContinue;
0856 
0857 #ifdef USE_CALLGRIND
0858   CALLGRIND_STOP_INSTRUMENTATION;
0859   CALLGRIND_DUMP_STATS;
0860 #endif
0861 }
0862 
0863 void MuScleFit::selectMuons(const edm::Event& event) {
0864   recMu1 = reco::Particle::LorentzVector(0, 0, 0, 0);
0865   recMu2 = reco::Particle::LorentzVector(0, 0, 0, 0);
0866 
0867   std::vector<MuScleFitMuon> muons;
0868   muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
0869   //  plotter->fillRec(muons); // @EM method already invoked inside MuScleFitMuonSelector::selectMuons()
0870 
0871   if (debug_ > 0) {
0872     std::cout << "[MuScleFit::selectMuons] Debugging muons collections after call to muonSelector_->selectMuons"
0873               << std::endl;
0874     int iMu = 0;
0875     for (std::vector<MuScleFitMuon>::const_iterator it = muons.begin(); it < muons.end(); ++it) {
0876       std::cout << "  - muon n. " << iMu << " = " << (*it) << std::endl;
0877       ++iMu;
0878     }
0879   }
0880 
0881   // Find the two muons from the resonance, and set ResFound bool
0882   // ------------------------------------------------------------
0883   std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes = MuScleFitUtils::findBestRecoRes(muons);
0884 
0885   if (MuScleFitUtils::ResFound) {
0886     if (debug_ > 0) {
0887       std::cout << std::setprecision(9) << "Pt after findbestrecores: " << (recMuFromBestRes.first).Pt() << " "
0888                 << (recMuFromBestRes.second).Pt() << std::endl;
0889       std::cout << "recMu1 = " << recMu1 << std::endl;
0890       std::cout << "recMu2 = " << recMu2 << std::endl;
0891     }
0892     recMu1 = recMuFromBestRes.first.p4();
0893     recMu2 = recMuFromBestRes.second.p4();
0894     recMuScleMu1 = recMuFromBestRes.first;
0895     recMuScleMu2 = recMuFromBestRes.second;
0896 
0897     if (debug_ > 0) {
0898       std::cout << "after recMu1 = " << recMu1 << std::endl;
0899       std::cout << "after recMu2 = " << recMu2 << std::endl;
0900       std::cout << "mu1.pt = " << recMu1.Pt() << std::endl;
0901       std::cout << "mu2.pt = " << recMu2.Pt() << std::endl;
0902       std::cout << "after recMuScleMu1 = " << recMuScleMu1 << std::endl;
0903       std::cout << "after recMuScleMu2 = " << recMuScleMu2 << std::endl;
0904     }
0905     MuScleFitUtils::SavedPair.push_back(std::make_pair(recMu1, recMu2));
0906     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(recMuScleMu1, recMuScleMu2));
0907   } else {
0908     MuScleFitUtils::SavedPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
0909     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(MuScleFitMuon(), MuScleFitMuon()));
0910   }
0911   // Save the events also in the external tree so that it can be saved late
0912 
0913   // Fetch extra information (per event)
0914   UInt_t the_NVtx(0);
0915   Int_t the_numPUvtx(0);
0916   Float_t the_TrueNumInteractions(0);
0917 
0918   // Fill pile-up related informations
0919   // --------------------------------
0920   edm::Handle<std::vector<PileupSummaryInfo> > puInfo;
0921   event.getByLabel(puInfoSrc_, puInfo);
0922   if (puInfo.isValid()) {
0923     std::vector<PileupSummaryInfo>::const_iterator PVI;
0924     for (PVI = puInfo->begin(); PVI != puInfo->end(); ++PVI) {
0925       int BX = PVI->getBunchCrossing();
0926       if (BX == 0) {  // "0" is the in-time crossing, negative values are the early crossings, positive are late
0927         the_TrueNumInteractions = PVI->getTrueNumInteractions();
0928         the_numPUvtx = PVI->getPU_NumInteractions();
0929       }
0930     }
0931   }
0932 
0933   edm::Handle<std::vector<reco::Vertex> > vertices;
0934   event.getByLabel(vertexSrc_, vertices);
0935   if (vertices.isValid()) {
0936     std::vector<reco::Vertex>::const_iterator itv;
0937     // now, count vertices
0938     for (itv = vertices->begin(); itv != vertices->end(); ++itv) {
0939       // require that the vertex meets certain criteria
0940       if (itv->ndof() < 5)
0941         continue;
0942       if (fabs(itv->z()) > 50.0)
0943         continue;
0944       if (fabs(itv->position().rho()) > 2.0)
0945         continue;
0946       ++the_NVtx;
0947     }
0948   }
0949 
0950   // get the MC event weight
0951   edm::Handle<GenEventInfoProduct> genEvtInfo;
0952   event.getByLabel("generator", genEvtInfo);
0953   double the_genEvtweight = 1.;
0954   if (genEvtInfo.isValid()) {
0955     the_genEvtweight = genEvtInfo->weight();
0956   }
0957 
0958   muonPairs_.push_back(MuonPair(
0959       MuScleFitUtils::SavedPairMuScleFitMuons.back().first,
0960       MuScleFitUtils::SavedPairMuScleFitMuons.back().second,
0961       MuScleFitEvent(
0962           event.run(), event.id().event(), the_genEvtweight, the_numPUvtx, the_TrueNumInteractions, the_NVtx)));
0963   // Fill the internal genPair tree from the external one
0964   if (MuScleFitUtils::speedup == false) {
0965     MuScleFitUtils::genPair.push_back(std::make_pair(genMuonPairs_.back().mu1.p4(), genMuonPairs_.back().mu2.p4()));
0966   }
0967 }
0968 
0969 void MuScleFit::selectMuons(const int maxEvents, const TString& treeFileName) {
0970   std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
0971   RootTreeHandler rootTreeHandler;
0972   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0973   if (MuScleFitUtils::speedup) {
0974     rootTreeHandler.readTree(
0975         maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPairMuScleFitMuons), theMuonType_, &evtRun);
0976   } else {
0977     rootTreeHandler.readTree(maxEvents,
0978                              inputRootTreeFileName_,
0979                              &(MuScleFitUtils::SavedPairMuScleFitMuons),
0980                              theMuonType_,
0981                              &evtRun,
0982                              &(MuScleFitUtils::genMuscleFitPair));
0983   }
0984   // Now loop on all the pairs and apply any smearing and bias if needed
0985   std::vector<std::pair<unsigned int, unsigned long long> >::iterator evtRunIt = evtRun.begin();
0986   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator it = MuScleFitUtils::SavedPairMuScleFitMuons.begin();
0987   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator genIt;
0988   if (MuScleFitUtils::speedup == false)
0989     genIt = MuScleFitUtils::genMuscleFitPair.begin();
0990   for (; it != MuScleFitUtils::SavedPairMuScleFitMuons.end(); ++it, ++evtRunIt) {
0991     // Apply any cut if requested
0992     // Note that cuts here are only applied to already selected muons. They should not be used unless
0993     // you are sure that the difference is negligible (e.g. the number of events with > 2 muons is negligible).
0994     double pt1 = it->first.pt();
0995     //std::cout << "pt1 = " << pt1 << std::endl;
0996     double pt2 = it->second.pt();
0997     //std::cout << "pt2 = " << pt2 << std::endl;
0998     double eta1 = it->first.eta();
0999     //std::cout << "eta1 = " << eta1 << std::endl;
1000     double eta2 = it->second.eta();
1001     //std::cout << "eta2 = " << eta2 << std::endl;
1002     // If they don't pass the cuts set to null vectors
1003     bool dontPass = false;
1004     bool eta1InFirstRange;
1005     bool eta2InFirstRange;
1006     bool eta1InSecondRange;
1007     bool eta2InSecondRange;
1008 
1009     int ch1 = it->first.charge();
1010     int ch2 = it->second.charge();
1011 
1012     if (MuScleFitUtils::separateRanges_) {
1013       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1014       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1015       eta1InSecondRange =
1016           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1017       eta2InSecondRange =
1018           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1019 
1020       // This is my logic, which should be erroneous, but certainly simpler...
1021       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1022             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1023             ((eta1InFirstRange && eta2InSecondRange && ch1 >= ch2) ||
1024              (eta1InSecondRange && eta2InFirstRange && ch1 < ch2))))
1025         dontPass = true;
1026     } else {
1027       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1028       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1029       eta1InSecondRange =
1030           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1031       eta2InSecondRange =
1032           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1033       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1034             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1035             (((eta1InFirstRange && !eta2InFirstRange) && (eta2InSecondRange && !eta1InSecondRange) && ch1 >= ch2) ||
1036              ((eta2InFirstRange && !eta1InFirstRange) && (eta1InSecondRange && !eta2InSecondRange) && ch1 < ch2))))
1037         dontPass = true;
1038     }
1039 
1040     // Additional check on deltaPhi
1041     double deltaPhi = MuScleFitUtils::deltaPhi(it->first.phi(), it->second.phi());
1042     if ((deltaPhi <= MuScleFitUtils::deltaPhiMinCut_) || (deltaPhi >= MuScleFitUtils::deltaPhiMaxCut_))
1043       dontPass = true;
1044 
1045     lorentzVector vec1 = it->first.p4();
1046     lorentzVector vec2 = it->second.p4();
1047     if (ch1 >= ch2) {
1048       lorentzVector vectemp = vec1;
1049       vec1 = vec2;
1050       vec2 = vectemp;
1051     }
1052 
1053     if (!dontPass) {
1054       // First is always mu-, second mu+
1055       if ((MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0)) {
1056         applySmearing(vec1);
1057         applyBias(vec1, -1);
1058         applySmearing(vec2);
1059         applyBias(vec2, 1);
1060       }
1061 
1062       MuScleFitUtils::SavedPair.push_back(std::make_pair(vec1, vec2));
1063     }
1064 
1065     //FIXME: we loose the additional information besides the 4-momenta
1066     muonPairs_.push_back(MuonPair(
1067         MuScleFitMuon(vec1, -1),
1068         MuScleFitMuon(vec2, +1),
1069         MuScleFitEvent(
1070             (*evtRunIt).first, (*evtRunIt).second, 0, 0, 0, 0))  // FIXME: order of event and run number mixed up!
1071     );
1072 
1073     // Fill the internal genPair tree from the external one
1074     if (!MuScleFitUtils::speedup) {
1075       MuScleFitUtils::genPair.push_back(std::make_pair(genIt->first.p4(), genIt->second.p4()));
1076       genMuonPairs_.push_back(GenMuonPair(genIt->first.p4(), genIt->second.p4(), 0));
1077       ++genIt;
1078     }
1079   }
1080   plotter->fillTreeRec(MuScleFitUtils::SavedPair);
1081   if (!(MuScleFitUtils::speedup)) {
1082     plotter->fillTreeGen(MuScleFitUtils::genPair);
1083   }
1084 }
1085 
1086 void MuScleFit::duringFastLoop() {
1087   // On loops>0 the two muons are directly obtained from the SavedMuon array
1088   // -----------------------------------------------------------------------
1089   MuScleFitUtils::ResFound = false;
1090   recMu1 = (MuScleFitUtils::SavedPair[iev].first);
1091   recMu2 = (MuScleFitUtils::SavedPair[iev].second);
1092 
1093   //std::cout << "iev = " << iev << ", recMu1 pt = " << recMu1.Pt() << ", recMu2 pt = " << recMu2.Pt() << std::endl;
1094 
1095   if (recMu1.Pt() > 0 && recMu2.Pt() > 0) {
1096     MuScleFitUtils::ResFound = true;
1097     if (debug_ > 0)
1098       std::cout << "Ev = " << iev << ": found muons in tree with Pt = " << recMu1.Pt() << " " << recMu2.Pt()
1099                 << std::endl;
1100   }
1101 
1102   if (debug_ > 0)
1103     std::cout << "About to start lik par correction and histo filling; ResFound is " << MuScleFitUtils::ResFound
1104               << std::endl;
1105   // If resonance found, do the hard work
1106   // ------------------------------------
1107   if (MuScleFitUtils::ResFound) {
1108     // Find weight and reference mass for this muon pair
1109     // -------------------------------------------------
1110     // The last parameter = true means that we want to use always the background window to compute the weight,
1111     // otherwise the probability will be filled only for the resonance region.
1112     double weight = MuScleFitUtils::computeWeight((recMu1 + recMu2).mass(), iev, true);
1113     if (debug_ > 0) {
1114       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": before correction     Pt1 = " << recMu1.Pt()
1115                 << " Pt2 = " << recMu2.Pt() << std::endl;
1116     }
1117     // For successive iterations, correct the muons only if the previous iteration was a scale fit.
1118     // --------------------------------------------------------------------------------------------
1119     if (loopCounter > 0) {
1120       if (MuScleFitUtils::doScaleFit[loopCounter - 1]) {
1121         recMu1 = (MuScleFitUtils::applyScale(recMu1, MuScleFitUtils::parvalue[loopCounter - 1], -1));
1122         recMu2 = (MuScleFitUtils::applyScale(recMu2, MuScleFitUtils::parvalue[loopCounter - 1], 1));
1123       }
1124     }
1125     if (debug_ > 0) {
1126       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": after correction      Pt1 = " << recMu1.Pt()
1127                 << " Pt2 = " << recMu2.Pt() << std::endl;
1128     }
1129 
1130     reco::Particle::LorentzVector bestRecRes(recMu1 + recMu2);
1131 
1132     //Fill histograms
1133     //------------------
1134 
1135     mapHisto_["hRecBestMu"]->Fill(recMu1, -1, weight);
1136     mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
1137     mapHisto_["hRecBestMu"]->Fill(recMu2, +1, weight);
1138     mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
1139     mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
1140     // Reconstructed resonance
1141     mapHisto_["hRecBestRes"]->Fill(bestRecRes, +1, weight);
1142     mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, +1, 1.);
1143     //     // Fill histogram of Res mass vs muon variables
1144     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu1, bestRecRes, -1);
1145     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu2, bestRecRes, +1);
1146     //     // Fill also the mass mu+/mu- comparisons
1147     //     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes);
1148 
1149     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, bestRecRes, -1, weight);
1150     mapHisto_["hRecBestResVSMu"]->Fill(recMu2, bestRecRes, +1, weight);
1151     // Fill also the mass mu+/mu- comparisons
1152     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes, weight);
1153 
1154     //-- rc 2010 filling histograms for mu+ /mu- ------
1155     //  mapHisto_["hRecBestResVSMuMinus"]->Fill (recMu1, bestRecRes, -1);
1156     // mapHisto_["hRecBestResVSMuPlus"]->Fill (recMu2, bestRecRes, +1);
1157 
1158     //-- rc 2010 filling histograms MassVsMuEtaPhi------
1159     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu1, bestRecRes,-1);
1160     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu2, bestRecRes,+1);
1161 
1162     // Fill histogram of Res mass vs Res variables
1163     // mapHisto_["hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
1164     mapHisto_["hRecBestResVSRes"]->Fill(bestRecRes, bestRecRes, +1, weight);
1165 
1166     std::vector<double>* parval;
1167     std::vector<double> initpar;
1168     // Store a pointer to the vector of parameters of the last iteration, or the initial
1169     // parameters if this is the first iteration
1170     if (loopCounter == 0) {
1171       initpar = MuScleFitUtils::parResol;
1172       initpar.insert(initpar.end(), MuScleFitUtils::parScale.begin(), MuScleFitUtils::parScale.end());
1173       initpar.insert(initpar.end(), MuScleFitUtils::parCrossSection.begin(), MuScleFitUtils::parCrossSection.end());
1174       initpar.insert(initpar.end(), MuScleFitUtils::parBgr.begin(), MuScleFitUtils::parBgr.end());
1175       parval = &initpar;
1176     } else {
1177       parval = &(MuScleFitUtils::parvalue[loopCounter - 1]);
1178     }
1179 
1180     //Compute pt resolution w.r.t generated and simulated muons
1181     //--------------------------------------------------------
1182     if (!MuScleFitUtils::speedup) {
1183       //first is always mu-, second is always mu+
1184       if (checkDeltaR(MuScleFitUtils::genPair[iev].first, recMu1)) {
1185         fillComparisonHistograms(MuScleFitUtils::genPair[iev].first, recMu1, "Gen", -1);
1186       }
1187       if (checkDeltaR(MuScleFitUtils::genPair[iev].second, recMu2)) {
1188         fillComparisonHistograms(MuScleFitUtils::genPair[iev].second, recMu2, "Gen", +1);
1189       }
1190       if (compareToSimTracks_) {
1191         //first is always mu-, second is always mu+
1192         if (checkDeltaR(MuScleFitUtils::simPair[iev].first, recMu1)) {
1193           fillComparisonHistograms(MuScleFitUtils::simPair[iev].first, recMu1, "Sim", -1);
1194         }
1195         if (checkDeltaR(MuScleFitUtils::simPair[iev].second, recMu2)) {
1196           fillComparisonHistograms(MuScleFitUtils::simPair[iev].second, recMu2, "Sim", +1);
1197         }
1198       }
1199     }
1200 
1201     // ATTENTION: this was done only when a matching was found. Moved it outside because, genInfo or not, we still want to see the resolution function
1202     // Fill also the resolution histogramsm using the resolution functions:
1203     // the parameters are those from the last iteration, as the muons up to this point have also the corrections from the same iteration.
1204     // Need to use a different array (ForVec), containing functors able to operate on std::vector<double>
1205     mapHisto_["hFunctionResolPt"]->Fill(
1206         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1207     mapHisto_["hFunctionResolCotgTheta"]->Fill(
1208         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1209     mapHisto_["hFunctionResolPhi"]->Fill(
1210         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1211     mapHisto_["hFunctionResolPt"]->Fill(
1212         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1213     mapHisto_["hFunctionResolCotgTheta"]->Fill(
1214         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1215     mapHisto_["hFunctionResolPhi"]->Fill(
1216         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1217 
1218     // Compute likelihood histograms
1219     // -----------------------------
1220     if (debug_ > 0)
1221       std::cout << "mass = " << bestRecRes.mass() << std::endl;
1222     if (weight != 0.) {
1223       double massResol;
1224       double prob;
1225       double deltalike;
1226       if (loopCounter == 0) {
1227         std::vector<double> initpar;
1228         initpar.reserve((int)(MuScleFitUtils::parResol.size()));
1229         for (int i = 0; i < (int)(MuScleFitUtils::parResol.size()); i++) {
1230           initpar.push_back(MuScleFitUtils::parResol[i]);
1231         }
1232         for (int i = 0; i < (int)(MuScleFitUtils::parScale.size()); i++) {
1233           initpar.push_back(MuScleFitUtils::parScale[i]);
1234         }
1235         //  for (int i=0; i<(int)(MuScleFitUtils::parCrossSection.size()); i++) {
1236         //    initpar.push_back(MuScleFitUtils::parCrossSection[i]);
1237         //  }
1238         MuScleFitUtils::crossSectionHandler->addParameters(initpar);
1239 
1240         for (int i = 0; i < (int)(MuScleFitUtils::parBgr.size()); i++) {
1241           initpar.push_back(MuScleFitUtils::parBgr[i]);
1242         }
1243         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, initpar);
1244         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(), massResol, initpar, true );
1245         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1246                                         bestRecRes.Eta(),
1247                                         bestRecRes.Rapidity(),
1248                                         massResol,
1249                                         initpar,
1250                                         true,
1251                                         recMu1.eta(),
1252                                         recMu2.eta());
1253       } else {
1254         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parvalue[loopCounter - 1]);
1255         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(),
1256         //                                       massResol, MuScleFitUtils::parvalue[loopCounter-1], true );
1257         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1258                                         bestRecRes.Eta(),
1259                                         bestRecRes.Rapidity(),
1260                                         massResol,
1261                                         MuScleFitUtils::parvalue[loopCounter - 1],
1262                                         true,
1263                                         recMu1.eta(),
1264                                         recMu2.eta());
1265       }
1266       if (debug_ > 0)
1267         std::cout << "inside weight: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1268       if (prob > 0) {
1269         if (debug_ > 0)
1270           std::cout << "inside prob: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1271 
1272         deltalike = log(prob) * weight;  // NB maximum likelihood --> deltalike is maximized
1273         mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike);
1274         mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike);
1275         mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike);
1276         mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike);
1277 
1278         double recoMass = (recMu1 + recMu2).mass();
1279         if (recoMass != 0) {
1280           // IMPORTANT: massResol is not a relative resolution
1281           mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1);
1282           mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1);
1283           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol / recoMass, -1);
1284           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol / recoMass, +1);
1285         }
1286 
1287         if (MuScleFitUtils::debugMassResol_) {
1288           mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1);
1289           mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1);
1290           mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1);
1291           mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1);
1292           mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1);
1293           mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1);
1294         }
1295 
1296         if (!MuScleFitUtils::speedup) {
1297           double genMass = (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second).mass();
1298           // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
1299           if (genMass != 0) {
1300             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].first),
1301                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1302                                            -1);
1303             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].second),
1304                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1305                                            +1);
1306             double diffMass = (recoMass - genMass) / genMass;
1307             // double diffMass = recoMass - genMass;
1308             // Fill if for both muons
1309             double pt1 = recMu1.pt();
1310             double eta1 = recMu1.eta();
1311             double pt2 = recMu2.pt();
1312             double eta2 = recMu2.eta();
1313             // This is to avoid nan
1314             if (diffMass == diffMass) {
1315               // Mass relative difference vs Pt and Eta. To be used to extract the true mass resolution
1316               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
1317               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
1318               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
1319               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
1320               // This is used for the covariance comparison
1321               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
1322               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
1323             } else {
1324               std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
1325             }
1326           }
1327           // Fill with mass resolution from resolution function
1328           double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
1329           mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
1330           mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
1331         }
1332 
1333         mapHisto_["hMass_P"]->Fill(bestRecRes.mass(), prob);
1334         if (debug_ > 0)
1335           std::cout << "mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1336         mapHisto_["hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
1337 
1338         mapHisto_["hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
1339         mapHisto_["hMassProbVsMu"]->Fill(recMu1, bestRecRes, -1, prob);
1340         mapHisto_["hMassProbVsMu"]->Fill(recMu2, bestRecRes, +1, prob);
1341         mapHisto_["hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
1342         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu1, bestRecRes, -1, prob);
1343         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu2, bestRecRes, +1, prob);
1344       }
1345     }
1346   }  // end if ResFound
1347 
1348   // Fill the pair
1349   // -------------
1350   if (loopCounter > 0) {
1351     if (debug_ > 0)
1352       std::cout << "[MuScleFit]: filling the pair" << std::endl;
1353     MuScleFitUtils::SavedPair[iev] = std::make_pair(recMu1, recMu2);
1354   }
1355 
1356   iev++;
1357   MuScleFitUtils::iev_++;
1358 
1359   // return kContinue;
1360 }
1361 
1362 bool MuScleFit::checkDeltaR(reco::Particle::LorentzVector& genMu, reco::Particle::LorentzVector& recMu) {
1363   //first is always mu-, second is always mu+
1364   double deltaR =
1365       sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
1366            ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
1367   if (deltaR < 0.01)
1368     return true;
1369   else if (debug_ > 0) {
1370     std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
1371               << " DOES NOT MATCH with generated muon from resonance: " << std::endl
1372               << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
1373   }
1374   return false;
1375 }
1376 
1377 void MuScleFit::fillComparisonHistograms(const reco::Particle::LorentzVector& genMu,
1378                                          const reco::Particle::LorentzVector& recMu,
1379                                          const std::string& inputName,
1380                                          const int charge) {
1381   std::string name(inputName + "VSMu");
1382   mapHisto_["hResolPt" + name]->Fill(recMu, (-genMu.Pt() + recMu.Pt()) / genMu.Pt(), charge);
1383   mapHisto_["hResolTheta" + name]->Fill(recMu, (-genMu.Theta() + recMu.Theta()), charge);
1384   mapHisto_["hResolCotgTheta" + name]->Fill(
1385       recMu, (-cos(genMu.Theta()) / sin(genMu.Theta()) + cos(recMu.Theta()) / sin(recMu.Theta())), charge);
1386   mapHisto_["hResolEta" + name]->Fill(recMu, (-genMu.Eta() + recMu.Eta()), charge);
1387   mapHisto_["hResolPhi" + name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);
1388 
1389   // Fill only if it was matched to a genMu and this muon is valid
1390   if ((genMu.Pt() != 0) && (recMu.Pt() != 0)) {
1391     mapHisto_["hPtRecoVsPt" + inputName]->Fill(genMu.Pt(), recMu.Pt());
1392   }
1393 }
1394 
1395 void MuScleFit::applySmearing(reco::Particle::LorentzVector& mu) {
1396   if (MuScleFitUtils::SmearType > 0) {
1397     mu = MuScleFitUtils::applySmearing(mu);
1398     if (debug_ > 0)
1399       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after smearing  Pt = " << mu.Pt() << std::endl;
1400   }
1401 }
1402 
1403 void MuScleFit::applyBias(reco::Particle::LorentzVector& mu, const int charge) {
1404   if (MuScleFitUtils::BiasType > 0) {
1405     mu = MuScleFitUtils::applyBias(mu, charge);
1406     if (debug_ > 0)
1407       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after bias      Pt = " << mu.Pt() << std::endl;
1408   }
1409 }
1410 
1411 // Simple method to check parameters consistency. It aborts the job if the parameters
1412 // are not consistent.
1413 void MuScleFit::checkParameters() {
1414   // Fits selection dimension check
1415   if (MuScleFitUtils::doResolFit.size() != maxLoopNumber) {
1416     std::cout << "[MuScleFit-Constructor]: wrong size of resolution fits selector = "
1417               << MuScleFitUtils::doResolFit.size() << std::endl;
1418     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1419     abort();
1420   }
1421   if (MuScleFitUtils::doScaleFit.size() != maxLoopNumber) {
1422     std::cout << "[MuScleFit-Constructor]: wrong size of scale fits selector = " << MuScleFitUtils::doScaleFit.size()
1423               << std::endl;
1424     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1425     abort();
1426   }
1427   if (MuScleFitUtils::doCrossSectionFit.size() != maxLoopNumber) {
1428     std::cout << "[MuScleFit-Constructor]: wrong size of cross section fits selector = "
1429               << MuScleFitUtils::doCrossSectionFit.size() << std::endl;
1430     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1431     abort();
1432   }
1433   if (MuScleFitUtils::doBackgroundFit.size() != maxLoopNumber) {
1434     std::cout << "[MuScleFit-Constructor]: wrong size of background fits selector = "
1435               << MuScleFitUtils::doBackgroundFit.size() << std::endl;
1436     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1437     abort();
1438   }
1439 
1440   // Bias parameters: dimension check
1441   // --------------------------------
1442   if ((MuScleFitUtils::BiasType == 1 && MuScleFitUtils::parBias.size() != 2) ||  // linear in pt
1443       (MuScleFitUtils::BiasType == 2 && MuScleFitUtils::parBias.size() != 2) ||  // linear in |eta|
1444       (MuScleFitUtils::BiasType == 3 && MuScleFitUtils::parBias.size() != 4) ||  // sinusoidal in phi
1445       (MuScleFitUtils::BiasType == 4 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and |eta|
1446       (MuScleFitUtils::BiasType == 5 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and sinusoidal in phi
1447       (MuScleFitUtils::BiasType == 6 &&
1448        MuScleFitUtils::parBias.size() != 3) ||  // linear in |eta| and sinusoidal in phi
1449       (MuScleFitUtils::BiasType == 7 &&
1450        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and |eta| and sinusoidal in phi
1451       (MuScleFitUtils::BiasType == 8 && MuScleFitUtils::parBias.size() != 4) ||   // linear in pt and parabolic in |eta|
1452       (MuScleFitUtils::BiasType == 9 && MuScleFitUtils::parBias.size() != 2) ||   // exponential in pt
1453       (MuScleFitUtils::BiasType == 10 && MuScleFitUtils::parBias.size() != 3) ||  // parabolic in pt
1454       (MuScleFitUtils::BiasType == 11 &&
1455        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and sin in phi with chg
1456       (MuScleFitUtils::BiasType == 12 &&
1457        MuScleFitUtils::parBias.size() != 6) ||  // linear in pt and para in plus sin in phi with chg
1458       (MuScleFitUtils::BiasType == 13 &&
1459        MuScleFitUtils::parBias.size() != 8) ||  // linear in pt and para in plus sin in phi with chg
1460       MuScleFitUtils::BiasType < 0 ||
1461       MuScleFitUtils::BiasType > 13) {
1462     std::cout << "[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
1463     abort();
1464   }
1465   // Smear parameters: dimension check
1466   // ---------------------------------
1467   if ((MuScleFitUtils::SmearType == 1 && MuScleFitUtils::parSmear.size() != 3) ||
1468       (MuScleFitUtils::SmearType == 2 && MuScleFitUtils::parSmear.size() != 4) ||
1469       (MuScleFitUtils::SmearType == 3 && MuScleFitUtils::parSmear.size() != 5) ||
1470       (MuScleFitUtils::SmearType == 4 && MuScleFitUtils::parSmear.size() != 6) ||
1471       (MuScleFitUtils::SmearType == 5 && MuScleFitUtils::parSmear.size() != 7) ||
1472       (MuScleFitUtils::SmearType == 6 && MuScleFitUtils::parSmear.size() != 16) ||
1473       (MuScleFitUtils::SmearType == 7 && !MuScleFitUtils::parSmear.empty()) || MuScleFitUtils::SmearType < 0 ||
1474       MuScleFitUtils::SmearType > 7) {
1475     std::cout << "[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
1476     abort();
1477   }
1478   // Protect against bad size of parameters
1479   // --------------------------------------
1480   if (MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolFix.size() ||
1481       MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolOrder.size()) {
1482     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
1483     abort();
1484   }
1485   if (MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleFix.size() ||
1486       MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleOrder.size()) {
1487     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
1488     abort();
1489   }
1490   if (MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionFix.size() ||
1491       MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionOrder.size()) {
1492     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1493     abort();
1494   }
1495   if (MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrFix.size() ||
1496       MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrOrder.size()) {
1497     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1498     abort();
1499   }
1500 
1501   // Protect against an incorrect number of resonances
1502   // -------------------------------------------------
1503   if (MuScleFitUtils::resfind.size() != 6) {
1504     std::cout << "[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
1505     abort();
1506   }
1507 }
1508 
1509 bool MuScleFit::selGlobalMuon(const pat::Muon* aMuon) {
1510   reco::TrackRef iTrack = aMuon->innerTrack();
1511   const reco::HitPattern& p = iTrack->hitPattern();
1512 
1513   reco::TrackRef gTrack = aMuon->globalTrack();
1514   const reco::HitPattern& q = gTrack->hitPattern();
1515 
1516   return (  //isMuonInAccept(aMuon) &&// no acceptance cuts!
1517       iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 && q.numberOfValidMuonHits() > 0 &&
1518       iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
1519       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
1520       fabs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
1521       fabs(iTrack->dz()) < 15.0);   //should be done w.r.t. PV!
1522 }
1523 
1524 bool MuScleFit::selTrackerMuon(const pat::Muon* aMuon) {
1525   reco::TrackRef iTrack = aMuon->innerTrack();
1526   const reco::HitPattern& p = iTrack->hitPattern();
1527 
1528   return (  //isMuonInAccept(aMuon) // no acceptance cuts!
1529       iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
1530       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
1531       fabs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
1532       fabs(iTrack->dz()) < 15.0);   //should be done w.r.t. PV!
1533 }
1534 
1535 #include "FWCore/Framework/interface/MakerMacros.h"
1536 DEFINE_FWK_LOOPER(MuScleFit);