Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-12 22:40:29

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   edm::EDGetTokenT<edm::TriggerResults> triggerResultsToken_;
0303   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0304   edm::EDGetTokenT<std::vector<PileupSummaryInfo> > puInfoToken_;
0305   edm::EDGetTokenT<GenEventInfoProduct> genEvtInfoToken_;
0306 
0307   std::string triggerResultsLabel_;
0308   std::string triggerResultsProcess_;
0309   std::vector<std::string> triggerPath_;
0310   bool negateTrigger_;
0311   bool saveAllToTree_;
0312 
0313   // input collections for PU related infos
0314   edm::InputTag puInfoSrc_;
0315   edm::InputTag vertexSrc_;
0316 
0317   std::unique_ptr<MuScleFitMuonSelector> muonSelector_;
0318 };
0319 
0320 template <typename T>
0321 std::vector<MuScleFitMuon> MuScleFit::fillMuonCollection(const std::vector<T>& tracks) {
0322   std::vector<MuScleFitMuon> muons;
0323   typename std::vector<T>::const_iterator track;
0324   for (track = tracks.begin(); track != tracks.end(); ++track) {
0325     reco::Particle::LorentzVector mu;
0326     mu = reco::Particle::LorentzVector(
0327         track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + MuScleFitUtils::mMu2));
0328     // Apply smearing if needed, and then bias
0329     // ---------------------------------------
0330     MuScleFitUtils::goodmuon++;
0331     if (debug_ > 0)
0332       std::cout << std::setprecision(9) << "Muon #" << MuScleFitUtils::goodmuon << ": initial value   Pt = " << mu.Pt()
0333                 << std::endl;
0334 
0335     applySmearing(mu);
0336     applyBias(mu, track->charge());
0337     if (debug_ > 0)
0338       std::cout << "track charge: " << track->charge() << std::endl;
0339 
0340     Double_t hitsTk = track->innerTrack()->hitPattern().numberOfValidTrackerHits();
0341     Double_t hitsMuon = track->innerTrack()->hitPattern().numberOfValidMuonHits();
0342     Double_t ptError = track->innerTrack()->ptError();
0343     MuScleFitMuon muon(mu, track->charge(), ptError, hitsTk, hitsMuon, false);
0344     if (debug_ > 0) {
0345       std::cout << "[MuScleFit::fillMuonCollection]" << std::endl;
0346       std::cout << "  muon = " << muon << std::endl;
0347     }
0348 
0349     // Store modified muon
0350     // -------------------
0351     muons.push_back(muon);
0352   }
0353   return muons;
0354 }
0355 
0356 template <typename T>
0357 void MuScleFit::takeSelectedMuonType(const T& muon, std::vector<reco::Track>& tracks) {
0358   // std::cout<<"muon "<<muon->isGlobalMuon()<<muon->isStandAloneMuon()<<muon->isTrackerMuon()<<std::endl;
0359   //NNBB: one muon can be of many kinds at once but with the theMuonType_ we are sure
0360   // to avoid double counting of the same muon
0361   if (muon->isGlobalMuon() && theMuonType_ == 1)
0362     tracks.push_back(*(muon->globalTrack()));
0363   else if (muon->isStandAloneMuon() && theMuonType_ == 2)
0364     tracks.push_back(*(muon->outerTrack()));
0365   else if (muon->isTrackerMuon() && theMuonType_ == 3)
0366     tracks.push_back(*(muon->innerTrack()));
0367 
0368   else if (theMuonType_ == 10 && !(muon->isStandAloneMuon()))  //particular case!!
0369     tracks.push_back(*(muon->innerTrack()));
0370   else if (theMuonType_ == 11 && muon->isGlobalMuon())
0371     tracks.push_back(*(muon->innerTrack()));
0372   else if (theMuonType_ == 13 && muon->isTrackerMuon())
0373     tracks.push_back(*(muon->innerTrack()));
0374 }
0375 
0376 // Constructor
0377 // -----------
0378 MuScleFit::MuScleFit(const edm::ParameterSet& pset) : MuScleFitBase(pset), totalEvents_(0) {
0379   MuScleFitUtils::debug = debug_;
0380   if (debug_ > 0)
0381     std::cout << "[MuScleFit]: Constructor" << std::endl;
0382 
0383   if ((theMuonType_ < -4 || theMuonType_ > 5) && theMuonType_ < 10) {
0384     std::cout << "[MuScleFit]: Unknown muon type! Aborting." << std::endl;
0385     abort();
0386   }
0387 
0388   loopCounter = 0;
0389 
0390   // Boundaries for h-function computation (to be improved!)
0391   // -------------------------------------------------------
0392   minResMass_hwindow[0] = 71.1876;  // 76.;
0393   maxResMass_hwindow[0] = 111.188;  // 106.;
0394   minResMass_hwindow[1] = 10.15;
0395   maxResMass_hwindow[1] = 10.55;
0396   minResMass_hwindow[2] = 9.8;
0397   maxResMass_hwindow[2] = 10.2;
0398   minResMass_hwindow[3] = 9.25;
0399   maxResMass_hwindow[3] = 9.65;
0400   minResMass_hwindow[4] = 3.58;
0401   maxResMass_hwindow[4] = 3.78;
0402   minResMass_hwindow[5] = 3.0;
0403   maxResMass_hwindow[5] = 3.2;
0404 
0405   // Max number of loops (if > 2 then try to minimize likelihood more than once)
0406   // ---------------------------------------------------------------------------
0407   maxLoopNumber = pset.getUntrackedParameter<int>("maxLoopNumber", 2);
0408   fastLoop = pset.getUntrackedParameter<bool>("FastLoop", true);
0409 
0410   // Selection of fits according to loop
0411   MuScleFitUtils::doResolFit = pset.getParameter<std::vector<int> >("doResolFit");
0412   MuScleFitUtils::doScaleFit = pset.getParameter<std::vector<int> >("doScaleFit");
0413   MuScleFitUtils::doCrossSectionFit = pset.getParameter<std::vector<int> >("doCrossSectionFit");
0414   MuScleFitUtils::doBackgroundFit = pset.getParameter<std::vector<int> >("doBackgroundFit");
0415 
0416   // Bias and smear types
0417   // --------------------
0418   int biasType = pset.getParameter<int>("BiasType");
0419   MuScleFitUtils::BiasType = biasType;
0420   // No error, the scale functions are used also for the bias
0421   MuScleFitUtils::biasFunction = scaleFunctionVecService(biasType);
0422   int smearType = pset.getParameter<int>("SmearType");
0423   MuScleFitUtils::SmearType = smearType;
0424   MuScleFitUtils::smearFunction = smearFunctionService(smearType);
0425 
0426   // Fit types
0427   // ---------
0428   int resolFitType = pset.getParameter<int>("ResolFitType");
0429   MuScleFitUtils::ResolFitType = resolFitType;
0430   MuScleFitUtils::resolutionFunction = resolutionFunctionService(resolFitType);
0431   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolFitType);
0432   int scaleType = pset.getParameter<int>("ScaleFitType");
0433   MuScleFitUtils::ScaleFitType = scaleType;
0434   MuScleFitUtils::scaleFunction = scaleFunctionService(scaleType);
0435   MuScleFitUtils::scaleFunctionForVec = scaleFunctionVecService(scaleType);
0436 
0437   // Initial parameters values
0438   // -------------------------
0439   MuScleFitUtils::parBias = pset.getParameter<std::vector<double> >("parBias");
0440   MuScleFitUtils::parSmear = pset.getParameter<std::vector<double> >("parSmear");
0441   MuScleFitUtils::parResol = pset.getParameter<std::vector<double> >("parResol");
0442   MuScleFitUtils::parResolStep =
0443       pset.getUntrackedParameter<std::vector<double> >("parResolStep", std::vector<double>());
0444   MuScleFitUtils::parResolMin = pset.getUntrackedParameter<std::vector<double> >("parResolMin", std::vector<double>());
0445   MuScleFitUtils::parResolMax = pset.getUntrackedParameter<std::vector<double> >("parResolMax", std::vector<double>());
0446   MuScleFitUtils::parScale = pset.getParameter<std::vector<double> >("parScale");
0447   MuScleFitUtils::parScaleStep =
0448       pset.getUntrackedParameter<std::vector<double> >("parScaleStep", std::vector<double>());
0449   MuScleFitUtils::parScaleMin = pset.getUntrackedParameter<std::vector<double> >("parScaleMin", std::vector<double>());
0450   MuScleFitUtils::parScaleMax = pset.getUntrackedParameter<std::vector<double> >("parScaleMax", std::vector<double>());
0451   MuScleFitUtils::parCrossSection = pset.getParameter<std::vector<double> >("parCrossSection");
0452   MuScleFitUtils::parBgr = pset.getParameter<std::vector<double> >("parBgr");
0453   MuScleFitUtils::parResolFix = pset.getParameter<std::vector<int> >("parResolFix");
0454   MuScleFitUtils::parScaleFix = pset.getParameter<std::vector<int> >("parScaleFix");
0455   MuScleFitUtils::parCrossSectionFix = pset.getParameter<std::vector<int> >("parCrossSectionFix");
0456   MuScleFitUtils::parBgrFix = pset.getParameter<std::vector<int> >("parBgrFix");
0457   MuScleFitUtils::parResolOrder = pset.getParameter<std::vector<int> >("parResolOrder");
0458   MuScleFitUtils::parScaleOrder = pset.getParameter<std::vector<int> >("parScaleOrder");
0459   MuScleFitUtils::parCrossSectionOrder = pset.getParameter<std::vector<int> >("parCrossSectionOrder");
0460   MuScleFitUtils::parBgrOrder = pset.getParameter<std::vector<int> >("parBgrOrder");
0461 
0462   MuScleFitUtils::resfind = pset.getParameter<std::vector<int> >("resfind");
0463   MuScleFitUtils::FitStrategy = pset.getParameter<int>("FitStrategy");
0464 
0465   // Option to skip unnecessary stuff
0466   // --------------------------------
0467   MuScleFitUtils::speedup = pset.getParameter<bool>("speedup");
0468 
0469   // Option to skip simTracks comparison
0470   compareToSimTracks_ = pset.getParameter<bool>("compareToSimTracks");
0471   simTracksCollection_ = pset.getUntrackedParameter<edm::InputTag>("SimTracksCollection", edm::InputTag("g4SimHits"));
0472 
0473   triggerResultsLabel_ = pset.getUntrackedParameter<std::string>("TriggerResultsLabel");
0474   triggerResultsProcess_ = pset.getUntrackedParameter<std::string>("TriggerResultsProcess");
0475 
0476   triggerResultsToken_ =
0477       consumes<edm::TriggerResults>(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()));
0478 
0479   triggerPath_ = pset.getUntrackedParameter<std::vector<std::string> >("TriggerPath");
0480   negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
0481   saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
0482 
0483   // input collections for PU related infos
0484   puInfoSrc_ = pset.getUntrackedParameter<edm::InputTag>("PileUpSummaryInfo");
0485   puInfoToken_ = consumes<std::vector<PileupSummaryInfo> >(puInfoSrc_);
0486 
0487   vertexSrc_ = pset.getUntrackedParameter<edm::InputTag>("PrimaryVertexCollection");
0488   vertexToken_ = consumes<reco::VertexCollection>(vertexSrc_);
0489 
0490   genEvtInfoToken_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0491 
0492   PATmuons_ = pset.getUntrackedParameter<bool>("PATmuons", false);
0493   genParticlesName_ = pset.getUntrackedParameter<std::string>("GenParticlesName", "genParticles");
0494 
0495   // Use the probability file or not. If not it will perform a simpler selection taking the muon pair with
0496   // invariant mass closer to the pdf value and will crash if some fit is attempted.
0497   MuScleFitUtils::useProbsFile_ = pset.getUntrackedParameter<bool>("UseProbsFile", true);
0498 
0499   // This must be set to true if using events generated with Sherpa
0500   MuScleFitUtils::sherpa_ = pset.getUntrackedParameter<bool>("Sherpa", false);
0501 
0502   MuScleFitUtils::rapidityBinsForZ_ = pset.getUntrackedParameter<bool>("RapidityBinsForZ", true);
0503 
0504   // Set the cuts on muons to be used in the fit
0505   MuScleFitUtils::separateRanges_ = pset.getUntrackedParameter<bool>("SeparateRanges", true);
0506   MuScleFitUtils::maxMuonPt_ = pset.getUntrackedParameter<double>("MaxMuonPt", 100000000.);
0507   MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
0508   MuScleFitUtils::minMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MinMuonEtaFirstRange", -6.);
0509   MuScleFitUtils::maxMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaFirstRange", 6.);
0510   MuScleFitUtils::minMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MinMuonEtaSecondRange", -100.);
0511   MuScleFitUtils::maxMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaSecondRange", 100.);
0512   MuScleFitUtils::deltaPhiMinCut_ = pset.getUntrackedParameter<double>("DeltaPhiMinCut", -100.);
0513   MuScleFitUtils::deltaPhiMaxCut_ = pset.getUntrackedParameter<double>("DeltaPhiMaxCut", 100.);
0514 
0515   MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("DebugMassResol", false);
0516   // MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents;
0517 
0518   // Check for parameters consistency
0519   // it will abort in case of errors.
0520   checkParameters();
0521 
0522   // Generate array of gaussian-distributed numbers for smearing
0523   // -----------------------------------------------------------
0524   if (MuScleFitUtils::SmearType > 0) {
0525     std::cout << "[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
0526     TF1 G("G", "[0]*exp(-0.5*pow(x,2))", -5., 5.);
0527     double norm = 1 / sqrt(2 * TMath::Pi());
0528     G.SetParameter(0, norm);
0529     for (int i = 0; i < 10000; i++) {
0530       for (int j = 0; j < 7; j++) {
0531         MuScleFitUtils::x[j][i] = G.GetRandom();
0532       }
0533     }
0534   }
0535   MuScleFitUtils::goodmuon = 0;
0536 
0537   if (theMuonType_ > 0 && theMuonType_ < 4) {
0538     MuScleFitUtils::MuonTypeForCheckMassWindow = theMuonType_ - 1;
0539     MuScleFitUtils::MuonType = theMuonType_ - 1;
0540   } else if (theMuonType_ == 0 || theMuonType_ == 4 || theMuonType_ == 5 || theMuonType_ >= 10 || theMuonType_ == -1 ||
0541              theMuonType_ == -2 || theMuonType_ == -3 || theMuonType_ == -4) {
0542     MuScleFitUtils::MuonTypeForCheckMassWindow = 2;
0543     MuScleFitUtils::MuonType = 2;
0544   } else {
0545     std::cout << "Wrong muon type " << theMuonType_ << std::endl;
0546     exit(1);
0547   }
0548 
0549   // When using standalone muons switch to the single Z pdf
0550   if (theMuonType_ == 2) {
0551     MuScleFitUtils::rapidityBinsForZ_ = false;
0552   }
0553 
0554   // Initialize ResMaxSigma And ResHalfWidth - 0 = global, 1 = SM, 2 = tracker
0555   // -------------------------------------------------------------------------
0556   MuScleFitUtils::massWindowHalfWidth[0][0] = 20.;
0557   MuScleFitUtils::massWindowHalfWidth[0][1] = 0.35;
0558   MuScleFitUtils::massWindowHalfWidth[0][2] = 0.35;
0559   MuScleFitUtils::massWindowHalfWidth[0][3] = 0.35;
0560   MuScleFitUtils::massWindowHalfWidth[0][4] = 0.2;
0561   MuScleFitUtils::massWindowHalfWidth[0][5] = 0.2;
0562   MuScleFitUtils::massWindowHalfWidth[1][0] = 20.;
0563   MuScleFitUtils::massWindowHalfWidth[1][1] = 0.35;
0564   MuScleFitUtils::massWindowHalfWidth[1][2] = 0.35;
0565   MuScleFitUtils::massWindowHalfWidth[1][3] = 0.35;
0566   MuScleFitUtils::massWindowHalfWidth[1][4] = 0.2;
0567   MuScleFitUtils::massWindowHalfWidth[1][5] = 0.2;
0568   MuScleFitUtils::massWindowHalfWidth[2][0] = 20.;
0569   MuScleFitUtils::massWindowHalfWidth[2][1] = 0.35;
0570   MuScleFitUtils::massWindowHalfWidth[2][2] = 0.35;
0571   MuScleFitUtils::massWindowHalfWidth[2][3] = 0.35;
0572   MuScleFitUtils::massWindowHalfWidth[2][4] = 0.2;
0573   MuScleFitUtils::massWindowHalfWidth[2][5] = 0.2;
0574 
0575   edm::ConsumesCollector iC = consumesCollector();
0576   muonSelector_ = std::make_unique<MuScleFitMuonSelector>(iC,
0577                                                           theMuonLabel_,
0578                                                           theMuonType_,
0579                                                           PATmuons_,
0580                                                           MuScleFitUtils::resfind,
0581                                                           MuScleFitUtils::speedup,
0582                                                           genParticlesName_,
0583                                                           compareToSimTracks_,
0584                                                           simTracksCollection_,
0585                                                           MuScleFitUtils::sherpa_,
0586                                                           debug_);
0587 
0588   MuScleFitUtils::backgroundHandler =
0589       new BackgroundHandler(pset.getParameter<std::vector<int> >("BgrFitType"),
0590                             pset.getParameter<std::vector<double> >("LeftWindowBorder"),
0591                             pset.getParameter<std::vector<double> >("RightWindowBorder"),
0592                             MuScleFitUtils::ResMass,
0593                             MuScleFitUtils::massWindowHalfWidth[MuScleFitUtils::MuonTypeForCheckMassWindow]);
0594 
0595   MuScleFitUtils::crossSectionHandler =
0596       new CrossSectionHandler(MuScleFitUtils::parCrossSection, MuScleFitUtils::resfind);
0597 
0598   // Build cross section scale factors
0599   // MuScleFitUtils::resfind
0600 
0601   MuScleFitUtils::normalizeLikelihoodByEventNumber_ =
0602       pset.getUntrackedParameter<bool>("NormalizeLikelihoodByEventNumber", true);
0603   if (debug_ > 0)
0604     std::cout << "End of MuScleFit constructor" << std::endl;
0605 
0606   inputRootTreeFileName_ = pset.getParameter<std::string>("InputRootTreeFileName");
0607   outputRootTreeFileName_ = pset.getParameter<std::string>("OutputRootTreeFileName");
0608   maxEventsFromRootTree_ = pset.getParameter<int>("MaxEventsFromRootTree");
0609 
0610   MuScleFitUtils::startWithSimplex_ = pset.getParameter<bool>("StartWithSimplex");
0611   MuScleFitUtils::computeMinosErrors_ = pset.getParameter<bool>("ComputeMinosErrors");
0612   MuScleFitUtils::minimumShapePlots_ = pset.getParameter<bool>("MinimumShapePlots");
0613 
0614   beginOfJobInConstructor();
0615 }
0616 
0617 // Destructor
0618 // ----------
0619 MuScleFit::~MuScleFit() {
0620   if (debug_ > 0)
0621     std::cout << "[MuScleFit]: Destructor" << std::endl;
0622   std::cout << "Total number of analyzed events = " << totalEvents_ << std::endl;
0623 
0624   if (!(outputRootTreeFileName_.empty())) {
0625     // 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_
0626     if (!(inputRootTreeFileName_.empty() && (int(MuScleFitUtils::SavedPair.size()) != totalEvents_))) {
0627       std::cout << "Saving muon pairs to root tree" << std::endl;
0628       RootTreeHandler rootTreeHandler;
0629       if (MuScleFitUtils::speedup) {
0630         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, 0, saveAllToTree_);
0631         if (debug_ > 0) {
0632           std::vector<MuonPair>::const_iterator it = muonPairs_.begin();
0633           std::cout << "[MuScleFit::~MuScleFit] (Destructor)" << std::endl;
0634           for (; it < muonPairs_.end(); ++it) {
0635             std::cout << "  Debugging pairs that are going to be written to file" << std::endl;
0636             std::cout << "  muon1 = " << it->mu1 << std::endl;
0637             std::cout << "  muon2 = " << it->mu2 << std::endl;
0638           }
0639         }
0640         rootTreeHandler.writeTree(outputRootTreeFileName_, &(muonPairs_), theMuonType_, nullptr, saveAllToTree_);
0641       } else {
0642         // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, &(MuScleFitUtils::genPair), saveAllToTree_ );
0643         rootTreeHandler.writeTree(
0644             outputRootTreeFileName_, &(muonPairs_), theMuonType_, &(genMuonPairs_), saveAllToTree_);
0645       }
0646     } else {
0647       std::cout << "ERROR: events in the vector = " << MuScleFitUtils::SavedPair.size()
0648                 << " != totalEvents = " << totalEvents_ << std::endl;
0649     }
0650   }
0651 }
0652 
0653 // Begin job
0654 // ---------
0655 void MuScleFit::beginOfJobInConstructor()
0656 // void MuScleFit::beginOfJob ()
0657 // void MuScleFit::beginOfJob( const edm::EventSetup& eventSetup )
0658 {
0659   if (debug_ > 0)
0660     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0661   //if(maxLoopNumber>1)
0662   if (MuScleFitUtils::useProbsFile_) {
0663     readProbabilityDistributionsFromFile();
0664   }
0665 
0666   if (debug_ > 0)
0667     std::cout << "[MuScleFit]: beginOfJob" << std::endl;
0668 
0669   // Create the root file
0670   // --------------------
0671   for (unsigned int i = 0; i < (maxLoopNumber); i++) {
0672     std::stringstream ss;
0673     ss << i;
0674     std::string rootFileName = ss.str() + "_" + theRootFileName_;
0675     if (theCompressionSettings_ > -1) {
0676       theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE", "", theCompressionSettings_));
0677     } else {
0678       theFiles_.push_back(new TFile(rootFileName.c_str(), "RECREATE"));
0679     }
0680   }
0681   if (debug_ > 0)
0682     std::cout << "[MuScleFit]: Root file created" << std::endl;
0683 
0684   std::cout << "creating plotter" << std::endl;
0685   plotter = new MuScleFitPlotter(theGenInfoRootFileName_);
0686   plotter->debug = debug_;
0687 }
0688 
0689 // End of job method
0690 // -----------------
0691 void MuScleFit::endOfJob() {
0692   if (debug_ > 0)
0693     std::cout << "[MuScleFit]: endOfJob" << std::endl;
0694 }
0695 
0696 // New loop
0697 // --------
0698 void MuScleFit::startingNewLoop(unsigned int iLoop) {
0699   if (debug_ > 0)
0700     std::cout << "[MuScleFit]: Starting loop # " << iLoop << std::endl;
0701 
0702   // Number of muons used
0703   // --------------------
0704   MuScleFitUtils::goodmuon = 0;
0705 
0706   // Counters for problem std::cout-ing
0707   // -----------------------------
0708   MuScleFitUtils::counter_resprob = 0;
0709 
0710   // Create the root file
0711   // --------------------
0712   fillHistoMap(theFiles_[iLoop], iLoop);
0713 
0714   loopCounter = iLoop;
0715   MuScleFitUtils::loopCounter = loopCounter;
0716 
0717   iev = 0;
0718   MuScleFitUtils::iev_ = 0;
0719 
0720   MuScleFitUtils::oldNormalization_ = 0;
0721 }
0722 
0723 // End of loop routine
0724 // -------------------
0725 edm::EDLooper::Status MuScleFit::endOfLoop(const edm::EventSetup& eventSetup, unsigned int iLoop) {
0726   unsigned int iFastLoop = 1;
0727 
0728   // Read the events from the root tree if requested
0729   if (!(inputRootTreeFileName_.empty())) {
0730     selectMuons(maxEventsFromRootTree_, inputRootTreeFileName_);
0731     // When reading from local file all the loops are done here
0732     totalEvents_ = MuScleFitUtils::SavedPair.size();
0733     iFastLoop = 0;
0734   } else {
0735     endOfFastLoop(iLoop);
0736   }
0737 
0738   // If a fastLoop is required we do all the remaining iterations here
0739   if (fastLoop == true) {
0740     for (; iFastLoop < maxLoopNumber; ++iFastLoop) {
0741       std::cout << "Starting fast loop number " << iFastLoop << std::endl;
0742 
0743       // In the first loop is called by the framework
0744       // if( iFastLoop > 0 ) {
0745       startingNewLoop(iFastLoop);
0746       // }
0747 
0748       // std::vector<std::pair<lorentzVector,lorentzVector> >::const_iterator it = MuScleFitUtils::SavedPair.begin();
0749       // for( ; it != SavedPair.end(); ++it ) {
0750       while (iev < totalEvents_) {
0751         if (iev % 50000 == 0) {
0752           std::cout << "Fast looping on event number " << iev << std::endl;
0753         }
0754         // This reads muons from SavedPair using iev to keep track of the event
0755         duringFastLoop();
0756       }
0757       std::cout << "End of fast loop number " << iFastLoop << ". Ran on " << iev << " events" << std::endl;
0758       endOfFastLoop(iFastLoop);
0759     }
0760   }
0761 
0762   if (iFastLoop >= maxLoopNumber - 1) {
0763     return kStop;
0764   } else {
0765     return kContinue;
0766   }
0767 }
0768 
0769 void MuScleFit::endOfFastLoop(const unsigned int iLoop) {
0770   // std::cout<< "Inside endOfFastLoop, iLoop = " << iLoop << " and loopCounter = " << loopCounter << std::endl;
0771 
0772   if (loopCounter == 0) {
0773     // plotter->writeHistoMap();
0774     // The destructor will call the writeHistoMap after the cd to the output file
0775     delete plotter;
0776   }
0777 
0778   std::cout << "Ending loop # " << iLoop << std::endl;
0779 
0780   // Write the histos to file
0781   // ------------------------
0782   // theFiles_[iLoop]->cd();
0783   writeHistoMap(iLoop);
0784 
0785   // Likelihood minimization to compute corrections
0786   // ----------------------------------------------
0787   // theFiles_[iLoop]->cd();
0788   TDirectory* likelihoodDir = theFiles_[iLoop]->mkdir("likelihood");
0789   likelihoodDir->cd();
0790   MuScleFitUtils::minimizeLikelihood();
0791 
0792   // ATTENTION, this was put BEFORE the minimizeLikelihood. Check for problems.
0793   theFiles_[iLoop]->Close();
0794   // ATTENTION: Check that this delete does not give any problem
0795   delete theFiles_[iLoop];
0796 
0797   // Clear the histos
0798   // ----------------
0799   clearHistoMap();
0800 }
0801 
0802 // Stuff to do during loop
0803 // -----------------------
0804 edm::EDLooper::Status MuScleFit::duringLoop(const edm::Event& event, const edm::EventSetup& eventSetup) {
0805   edm::Handle<edm::TriggerResults> triggerResults = event.getHandle(triggerResultsToken_);
0806   bool isFired = false;
0807 
0808   if (triggerPath_[0].empty())
0809     isFired = true;
0810   else if (triggerPath_[0] == "All") {
0811     isFired = triggerResults->accept();
0812     if (debug_ > 0)
0813       std::cout << "Trigger " << isFired << std::endl;
0814   } else {
0815     bool changed;
0816     HLTConfigProvider hltConfig;
0817     hltConfig.init(event.getRun(), eventSetup, triggerResultsProcess_, changed);
0818 
0819     const edm::TriggerNames& triggerNames = event.triggerNames(*triggerResults);
0820 
0821     for (unsigned i = 0; i < triggerNames.size(); i++) {
0822       const std::string& hltName = triggerNames.triggerName(i);
0823 
0824       // match the path in the pset with the true name of the trigger
0825       for (unsigned int ipath = 0; ipath < triggerPath_.size(); ipath++) {
0826         if (hltName.find(triggerPath_[ipath]) != std::string::npos) {
0827           unsigned int triggerIndex(hltConfig.triggerIndex(hltName));
0828 
0829           // triggerIndex must be less than the size of HLTR or you get a CMSException: _M_range_check
0830           if (triggerIndex < triggerResults->size()) {
0831             isFired = triggerResults->accept(triggerIndex);
0832             if (debug_ > 0)
0833               std::cout << triggerPath_[ipath] << " " << hltName << " " << isFired << std::endl;
0834           }
0835         }  // end if (matching the path in the pset with the true trigger name
0836       }
0837     }
0838   }
0839 
0840   if (negateTrigger_ && isFired)
0841     return kContinue;
0842   else if (!(negateTrigger_) && !isFired)
0843     return kContinue;
0844 
0845 #ifdef USE_CALLGRIND
0846   CALLGRIND_START_INSTRUMENTATION;
0847 #endif
0848 
0849   if (debug_ > 0) {
0850     std::cout << "[MuScleFit-duringLoop]: loopCounter = " << loopCounter << " Run: " << event.id().run()
0851               << " Event: " << event.id().event() << std::endl;
0852   }
0853 
0854   // On the first iteration we read the bank, otherwise we fetch the information from the muon tree
0855   // ------------------------------------ Important Note --------------------------------------- //
0856   // The fillMuonCollection method applies any smearing or bias to the muons, so we NEVER use
0857   // unbiased muons.
0858   // ----------------------------------------------------------------------------------------------
0859   if (loopCounter == 0) {
0860     if (!fastLoop || inputRootTreeFileName_.empty()) {
0861       if (debug_ > 0)
0862         std::cout << "Reading from edm event" << std::endl;
0863       selectMuons(event);
0864       duringFastLoop();
0865       ++totalEvents_;
0866     }
0867   }
0868 
0869   return kContinue;
0870 
0871 #ifdef USE_CALLGRIND
0872   CALLGRIND_STOP_INSTRUMENTATION;
0873   CALLGRIND_DUMP_STATS;
0874 #endif
0875 }
0876 
0877 void MuScleFit::selectMuons(const edm::Event& event) {
0878   recMu1 = reco::Particle::LorentzVector(0, 0, 0, 0);
0879   recMu2 = reco::Particle::LorentzVector(0, 0, 0, 0);
0880 
0881   std::vector<MuScleFitMuon> muons;
0882   muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
0883   //  plotter->fillRec(muons); // @EM method already invoked inside MuScleFitMuonSelector::selectMuons()
0884 
0885   if (debug_ > 0) {
0886     std::cout << "[MuScleFit::selectMuons] Debugging muons collections after call to muonSelector_->selectMuons"
0887               << std::endl;
0888     int iMu = 0;
0889     for (std::vector<MuScleFitMuon>::const_iterator it = muons.begin(); it < muons.end(); ++it) {
0890       std::cout << "  - muon n. " << iMu << " = " << (*it) << std::endl;
0891       ++iMu;
0892     }
0893   }
0894 
0895   // Find the two muons from the resonance, and set ResFound bool
0896   // ------------------------------------------------------------
0897   std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes = MuScleFitUtils::findBestRecoRes(muons);
0898 
0899   if (MuScleFitUtils::ResFound) {
0900     if (debug_ > 0) {
0901       std::cout << std::setprecision(9) << "Pt after findbestrecores: " << (recMuFromBestRes.first).Pt() << " "
0902                 << (recMuFromBestRes.second).Pt() << std::endl;
0903       std::cout << "recMu1 = " << recMu1 << std::endl;
0904       std::cout << "recMu2 = " << recMu2 << std::endl;
0905     }
0906     recMu1 = recMuFromBestRes.first.p4();
0907     recMu2 = recMuFromBestRes.second.p4();
0908     recMuScleMu1 = recMuFromBestRes.first;
0909     recMuScleMu2 = recMuFromBestRes.second;
0910 
0911     if (debug_ > 0) {
0912       std::cout << "after recMu1 = " << recMu1 << std::endl;
0913       std::cout << "after recMu2 = " << recMu2 << std::endl;
0914       std::cout << "mu1.pt = " << recMu1.Pt() << std::endl;
0915       std::cout << "mu2.pt = " << recMu2.Pt() << std::endl;
0916       std::cout << "after recMuScleMu1 = " << recMuScleMu1 << std::endl;
0917       std::cout << "after recMuScleMu2 = " << recMuScleMu2 << std::endl;
0918     }
0919     MuScleFitUtils::SavedPair.push_back(std::make_pair(recMu1, recMu2));
0920     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(recMuScleMu1, recMuScleMu2));
0921   } else {
0922     MuScleFitUtils::SavedPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
0923     MuScleFitUtils::SavedPairMuScleFitMuons.push_back(std::make_pair(MuScleFitMuon(), MuScleFitMuon()));
0924   }
0925   // Save the events also in the external tree so that it can be saved late
0926 
0927   // Fetch extra information (per event)
0928   UInt_t the_NVtx(0);
0929   Int_t the_numPUvtx(0);
0930   Float_t the_TrueNumInteractions(0);
0931 
0932   // Fill pile-up related informations
0933   // --------------------------------
0934   edm::Handle<std::vector<PileupSummaryInfo> > puInfo = event.getHandle(puInfoToken_);
0935   if (puInfo.isValid()) {
0936     std::vector<PileupSummaryInfo>::const_iterator PVI;
0937     for (PVI = puInfo->begin(); PVI != puInfo->end(); ++PVI) {
0938       int BX = PVI->getBunchCrossing();
0939       if (BX == 0) {  // "0" is the in-time crossing, negative values are the early crossings, positive are late
0940         the_TrueNumInteractions = PVI->getTrueNumInteractions();
0941         the_numPUvtx = PVI->getPU_NumInteractions();
0942       }
0943     }
0944   }
0945 
0946   edm::Handle<std::vector<reco::Vertex> > vertices = event.getHandle(vertexToken_);
0947   if (vertices.isValid()) {
0948     std::vector<reco::Vertex>::const_iterator itv;
0949     // now, count vertices
0950     for (itv = vertices->begin(); itv != vertices->end(); ++itv) {
0951       // require that the vertex meets certain criteria
0952       if (itv->ndof() < 5)
0953         continue;
0954       if (fabs(itv->z()) > 50.0)
0955         continue;
0956       if (fabs(itv->position().rho()) > 2.0)
0957         continue;
0958       ++the_NVtx;
0959     }
0960   }
0961 
0962   // get the MC event weight
0963   edm::Handle<GenEventInfoProduct> genEvtInfo = event.getHandle(genEvtInfoToken_);
0964   double the_genEvtweight = 1.;
0965   if (genEvtInfo.isValid()) {
0966     the_genEvtweight = genEvtInfo->weight();
0967   }
0968 
0969   muonPairs_.push_back(MuonPair(
0970       MuScleFitUtils::SavedPairMuScleFitMuons.back().first,
0971       MuScleFitUtils::SavedPairMuScleFitMuons.back().second,
0972       MuScleFitEvent(
0973           event.run(), event.id().event(), the_genEvtweight, the_numPUvtx, the_TrueNumInteractions, the_NVtx)));
0974   // Fill the internal genPair tree from the external one
0975   if (MuScleFitUtils::speedup == false) {
0976     MuScleFitUtils::genPair.push_back(std::make_pair(genMuonPairs_.back().mu1.p4(), genMuonPairs_.back().mu2.p4()));
0977   }
0978 }
0979 
0980 void MuScleFit::selectMuons(const int maxEvents, const TString& treeFileName) {
0981   std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
0982   RootTreeHandler rootTreeHandler;
0983   std::vector<std::pair<unsigned int, unsigned long long> > evtRun;
0984   if (MuScleFitUtils::speedup) {
0985     rootTreeHandler.readTree(
0986         maxEvents, inputRootTreeFileName_, &(MuScleFitUtils::SavedPairMuScleFitMuons), theMuonType_, &evtRun);
0987   } else {
0988     rootTreeHandler.readTree(maxEvents,
0989                              inputRootTreeFileName_,
0990                              &(MuScleFitUtils::SavedPairMuScleFitMuons),
0991                              theMuonType_,
0992                              &evtRun,
0993                              &(MuScleFitUtils::genMuscleFitPair));
0994   }
0995   // Now loop on all the pairs and apply any smearing and bias if needed
0996   std::vector<std::pair<unsigned int, unsigned long long> >::iterator evtRunIt = evtRun.begin();
0997   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator it = MuScleFitUtils::SavedPairMuScleFitMuons.begin();
0998   std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >::iterator genIt;
0999   if (MuScleFitUtils::speedup == false)
1000     genIt = MuScleFitUtils::genMuscleFitPair.begin();
1001   for (; it != MuScleFitUtils::SavedPairMuScleFitMuons.end(); ++it, ++evtRunIt) {
1002     // Apply any cut if requested
1003     // Note that cuts here are only applied to already selected muons. They should not be used unless
1004     // you are sure that the difference is negligible (e.g. the number of events with > 2 muons is negligible).
1005     double pt1 = it->first.pt();
1006     //std::cout << "pt1 = " << pt1 << std::endl;
1007     double pt2 = it->second.pt();
1008     //std::cout << "pt2 = " << pt2 << std::endl;
1009     double eta1 = it->first.eta();
1010     //std::cout << "eta1 = " << eta1 << std::endl;
1011     double eta2 = it->second.eta();
1012     //std::cout << "eta2 = " << eta2 << std::endl;
1013     // If they don't pass the cuts set to null vectors
1014     bool dontPass = false;
1015     bool eta1InFirstRange;
1016     bool eta2InFirstRange;
1017     bool eta1InSecondRange;
1018     bool eta2InSecondRange;
1019 
1020     int ch1 = it->first.charge();
1021     int ch2 = it->second.charge();
1022 
1023     if (MuScleFitUtils::separateRanges_) {
1024       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1025       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1026       eta1InSecondRange =
1027           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1028       eta2InSecondRange =
1029           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1030 
1031       // This is my logic, which should be erroneous, but certainly simpler...
1032       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1033             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1034             ((eta1InFirstRange && eta2InSecondRange && ch1 >= ch2) ||
1035              (eta1InSecondRange && eta2InFirstRange && ch1 < ch2))))
1036         dontPass = true;
1037     } else {
1038       eta1InFirstRange = eta1 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta1 < MuScleFitUtils::maxMuonEtaFirstRange_;
1039       eta2InFirstRange = eta2 >= MuScleFitUtils::minMuonEtaFirstRange_ && eta2 < MuScleFitUtils::maxMuonEtaFirstRange_;
1040       eta1InSecondRange =
1041           eta1 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta1 < MuScleFitUtils::maxMuonEtaSecondRange_;
1042       eta2InSecondRange =
1043           eta2 >= MuScleFitUtils::minMuonEtaSecondRange_ && eta2 < MuScleFitUtils::maxMuonEtaSecondRange_;
1044       if (!(pt1 >= MuScleFitUtils::minMuonPt_ && pt1 < MuScleFitUtils::maxMuonPt_ &&
1045             pt2 >= MuScleFitUtils::minMuonPt_ && pt2 < MuScleFitUtils::maxMuonPt_ &&
1046             (((eta1InFirstRange && !eta2InFirstRange) && (eta2InSecondRange && !eta1InSecondRange) && ch1 >= ch2) ||
1047              ((eta2InFirstRange && !eta1InFirstRange) && (eta1InSecondRange && !eta2InSecondRange) && ch1 < ch2))))
1048         dontPass = true;
1049     }
1050 
1051     // Additional check on deltaPhi
1052     double deltaPhi = MuScleFitUtils::deltaPhi(it->first.phi(), it->second.phi());
1053     if ((deltaPhi <= MuScleFitUtils::deltaPhiMinCut_) || (deltaPhi >= MuScleFitUtils::deltaPhiMaxCut_))
1054       dontPass = true;
1055 
1056     lorentzVector vec1 = it->first.p4();
1057     lorentzVector vec2 = it->second.p4();
1058     if (ch1 >= ch2) {
1059       lorentzVector vectemp = vec1;
1060       vec1 = vec2;
1061       vec2 = vectemp;
1062     }
1063 
1064     if (!dontPass) {
1065       // First is always mu-, second mu+
1066       if ((MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0)) {
1067         applySmearing(vec1);
1068         applyBias(vec1, -1);
1069         applySmearing(vec2);
1070         applyBias(vec2, 1);
1071       }
1072 
1073       MuScleFitUtils::SavedPair.push_back(std::make_pair(vec1, vec2));
1074     }
1075 
1076     //FIXME: we loose the additional information besides the 4-momenta
1077     muonPairs_.push_back(MuonPair(
1078         MuScleFitMuon(vec1, -1),
1079         MuScleFitMuon(vec2, +1),
1080         MuScleFitEvent(
1081             (*evtRunIt).first, (*evtRunIt).second, 0, 0, 0, 0))  // FIXME: order of event and run number mixed up!
1082     );
1083 
1084     // Fill the internal genPair tree from the external one
1085     if (!MuScleFitUtils::speedup) {
1086       MuScleFitUtils::genPair.push_back(std::make_pair(genIt->first.p4(), genIt->second.p4()));
1087       genMuonPairs_.push_back(GenMuonPair(genIt->first.p4(), genIt->second.p4(), 0));
1088       ++genIt;
1089     }
1090   }
1091   plotter->fillTreeRec(MuScleFitUtils::SavedPair);
1092   if (!(MuScleFitUtils::speedup)) {
1093     plotter->fillTreeGen(MuScleFitUtils::genPair);
1094   }
1095 }
1096 
1097 void MuScleFit::duringFastLoop() {
1098   // On loops>0 the two muons are directly obtained from the SavedMuon array
1099   // -----------------------------------------------------------------------
1100   MuScleFitUtils::ResFound = false;
1101   recMu1 = (MuScleFitUtils::SavedPair[iev].first);
1102   recMu2 = (MuScleFitUtils::SavedPair[iev].second);
1103 
1104   //std::cout << "iev = " << iev << ", recMu1 pt = " << recMu1.Pt() << ", recMu2 pt = " << recMu2.Pt() << std::endl;
1105 
1106   if (recMu1.Pt() > 0 && recMu2.Pt() > 0) {
1107     MuScleFitUtils::ResFound = true;
1108     if (debug_ > 0)
1109       std::cout << "Ev = " << iev << ": found muons in tree with Pt = " << recMu1.Pt() << " " << recMu2.Pt()
1110                 << std::endl;
1111   }
1112 
1113   if (debug_ > 0)
1114     std::cout << "About to start lik par correction and histo filling; ResFound is " << MuScleFitUtils::ResFound
1115               << std::endl;
1116   // If resonance found, do the hard work
1117   // ------------------------------------
1118   if (MuScleFitUtils::ResFound) {
1119     // Find weight and reference mass for this muon pair
1120     // -------------------------------------------------
1121     // The last parameter = true means that we want to use always the background window to compute the weight,
1122     // otherwise the probability will be filled only for the resonance region.
1123     double weight = MuScleFitUtils::computeWeight((recMu1 + recMu2).mass(), iev, true);
1124     if (debug_ > 0) {
1125       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": before correction     Pt1 = " << recMu1.Pt()
1126                 << " Pt2 = " << recMu2.Pt() << std::endl;
1127     }
1128     // For successive iterations, correct the muons only if the previous iteration was a scale fit.
1129     // --------------------------------------------------------------------------------------------
1130     if (loopCounter > 0) {
1131       if (MuScleFitUtils::doScaleFit[loopCounter - 1]) {
1132         recMu1 = (MuScleFitUtils::applyScale(recMu1, MuScleFitUtils::parvalue[loopCounter - 1], -1));
1133         recMu2 = (MuScleFitUtils::applyScale(recMu2, MuScleFitUtils::parvalue[loopCounter - 1], 1));
1134       }
1135     }
1136     if (debug_ > 0) {
1137       std::cout << "Loop #" << loopCounter << "Event #" << iev << ": after correction      Pt1 = " << recMu1.Pt()
1138                 << " Pt2 = " << recMu2.Pt() << std::endl;
1139     }
1140 
1141     reco::Particle::LorentzVector bestRecRes(recMu1 + recMu2);
1142 
1143     //Fill histograms
1144     //------------------
1145 
1146     mapHisto_["hRecBestMu"]->Fill(recMu1, -1, weight);
1147     mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
1148     mapHisto_["hRecBestMu"]->Fill(recMu2, +1, weight);
1149     mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
1150     mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
1151     // Reconstructed resonance
1152     mapHisto_["hRecBestRes"]->Fill(bestRecRes, +1, weight);
1153     mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, +1, 1.);
1154     //     // Fill histogram of Res mass vs muon variables
1155     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu1, bestRecRes, -1);
1156     //     mapHisto_["hRecBestResVSMu"]->Fill (recMu2, bestRecRes, +1);
1157     //     // Fill also the mass mu+/mu- comparisons
1158     //     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes);
1159 
1160     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, bestRecRes, -1, weight);
1161     mapHisto_["hRecBestResVSMu"]->Fill(recMu2, bestRecRes, +1, weight);
1162     // Fill also the mass mu+/mu- comparisons
1163     mapHisto_["hRecBestResVSMu"]->Fill(recMu1, recMu2, bestRecRes, weight);
1164 
1165     //-- rc 2010 filling histograms for mu+ /mu- ------
1166     //  mapHisto_["hRecBestResVSMuMinus"]->Fill (recMu1, bestRecRes, -1);
1167     // mapHisto_["hRecBestResVSMuPlus"]->Fill (recMu2, bestRecRes, +1);
1168 
1169     //-- rc 2010 filling histograms MassVsMuEtaPhi------
1170     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu1, bestRecRes,-1);
1171     //  mapHisto_["hRecBestResVSMuEtaPhi"]->Fill (recMu2, bestRecRes,+1);
1172 
1173     // Fill histogram of Res mass vs Res variables
1174     // mapHisto_["hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
1175     mapHisto_["hRecBestResVSRes"]->Fill(bestRecRes, bestRecRes, +1, weight);
1176 
1177     std::vector<double>* parval;
1178     std::vector<double> initpar;
1179     // Store a pointer to the vector of parameters of the last iteration, or the initial
1180     // parameters if this is the first iteration
1181     if (loopCounter == 0) {
1182       initpar = MuScleFitUtils::parResol;
1183       initpar.insert(initpar.end(), MuScleFitUtils::parScale.begin(), MuScleFitUtils::parScale.end());
1184       initpar.insert(initpar.end(), MuScleFitUtils::parCrossSection.begin(), MuScleFitUtils::parCrossSection.end());
1185       initpar.insert(initpar.end(), MuScleFitUtils::parBgr.begin(), MuScleFitUtils::parBgr.end());
1186       parval = &initpar;
1187     } else {
1188       parval = &(MuScleFitUtils::parvalue[loopCounter - 1]);
1189     }
1190 
1191     //Compute pt resolution w.r.t generated and simulated muons
1192     //--------------------------------------------------------
1193     if (!MuScleFitUtils::speedup) {
1194       //first is always mu-, second is always mu+
1195       if (checkDeltaR(MuScleFitUtils::genPair[iev].first, recMu1)) {
1196         fillComparisonHistograms(MuScleFitUtils::genPair[iev].first, recMu1, "Gen", -1);
1197       }
1198       if (checkDeltaR(MuScleFitUtils::genPair[iev].second, recMu2)) {
1199         fillComparisonHistograms(MuScleFitUtils::genPair[iev].second, recMu2, "Gen", +1);
1200       }
1201       if (compareToSimTracks_) {
1202         //first is always mu-, second is always mu+
1203         if (checkDeltaR(MuScleFitUtils::simPair[iev].first, recMu1)) {
1204           fillComparisonHistograms(MuScleFitUtils::simPair[iev].first, recMu1, "Sim", -1);
1205         }
1206         if (checkDeltaR(MuScleFitUtils::simPair[iev].second, recMu2)) {
1207           fillComparisonHistograms(MuScleFitUtils::simPair[iev].second, recMu2, "Sim", +1);
1208         }
1209       }
1210     }
1211 
1212     // 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
1213     // Fill also the resolution histogramsm using the resolution functions:
1214     // the parameters are those from the last iteration, as the muons up to this point have also the corrections from the same iteration.
1215     // Need to use a different array (ForVec), containing functors able to operate on std::vector<double>
1216     mapHisto_["hFunctionResolPt"]->Fill(
1217         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1218     mapHisto_["hFunctionResolCotgTheta"]->Fill(
1219         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1220     mapHisto_["hFunctionResolPhi"]->Fill(
1221         recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), *parval), -1);
1222     mapHisto_["hFunctionResolPt"]->Fill(
1223         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1224     mapHisto_["hFunctionResolCotgTheta"]->Fill(
1225         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1226     mapHisto_["hFunctionResolPhi"]->Fill(
1227         recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), *parval), +1);
1228 
1229     // Compute likelihood histograms
1230     // -----------------------------
1231     if (debug_ > 0)
1232       std::cout << "mass = " << bestRecRes.mass() << std::endl;
1233     if (weight != 0.) {
1234       double massResol;
1235       double prob;
1236       double deltalike;
1237       if (loopCounter == 0) {
1238         std::vector<double> initpar;
1239         initpar.reserve((int)(MuScleFitUtils::parResol.size()));
1240         for (int i = 0; i < (int)(MuScleFitUtils::parResol.size()); i++) {
1241           initpar.push_back(MuScleFitUtils::parResol[i]);
1242         }
1243         for (int i = 0; i < (int)(MuScleFitUtils::parScale.size()); i++) {
1244           initpar.push_back(MuScleFitUtils::parScale[i]);
1245         }
1246         //  for (int i=0; i<(int)(MuScleFitUtils::parCrossSection.size()); i++) {
1247         //    initpar.push_back(MuScleFitUtils::parCrossSection[i]);
1248         //  }
1249         MuScleFitUtils::crossSectionHandler->addParameters(initpar);
1250 
1251         for (int i = 0; i < (int)(MuScleFitUtils::parBgr.size()); i++) {
1252           initpar.push_back(MuScleFitUtils::parBgr[i]);
1253         }
1254         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, initpar);
1255         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(), massResol, initpar, true );
1256         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1257                                         bestRecRes.Eta(),
1258                                         bestRecRes.Rapidity(),
1259                                         massResol,
1260                                         initpar,
1261                                         true,
1262                                         recMu1.eta(),
1263                                         recMu2.eta());
1264       } else {
1265         massResol = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parvalue[loopCounter - 1]);
1266         // prob      = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(),
1267         //                                       massResol, MuScleFitUtils::parvalue[loopCounter-1], true );
1268         prob = MuScleFitUtils::massProb(bestRecRes.mass(),
1269                                         bestRecRes.Eta(),
1270                                         bestRecRes.Rapidity(),
1271                                         massResol,
1272                                         MuScleFitUtils::parvalue[loopCounter - 1],
1273                                         true,
1274                                         recMu1.eta(),
1275                                         recMu2.eta());
1276       }
1277       if (debug_ > 0)
1278         std::cout << "inside weight: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1279       if (prob > 0) {
1280         if (debug_ > 0)
1281           std::cout << "inside prob: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1282 
1283         deltalike = log(prob) * weight;  // NB maximum likelihood --> deltalike is maximized
1284         mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike);
1285         mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike);
1286         mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike);
1287         mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike);
1288 
1289         double recoMass = (recMu1 + recMu2).mass();
1290         if (recoMass != 0) {
1291           // IMPORTANT: massResol is not a relative resolution
1292           mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1);
1293           mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1);
1294           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol / recoMass, -1);
1295           mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol / recoMass, +1);
1296         }
1297 
1298         if (MuScleFitUtils::debugMassResol_) {
1299           mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1);
1300           mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1);
1301           mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1);
1302           mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1);
1303           mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1);
1304           mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1);
1305         }
1306 
1307         if (!MuScleFitUtils::speedup) {
1308           double genMass = (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second).mass();
1309           // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
1310           if (genMass != 0) {
1311             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].first),
1312                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1313                                            -1);
1314             mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].second),
1315                                            (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second),
1316                                            +1);
1317             double diffMass = (recoMass - genMass) / genMass;
1318             // double diffMass = recoMass - genMass;
1319             // Fill if for both muons
1320             double pt1 = recMu1.pt();
1321             double eta1 = recMu1.eta();
1322             double pt2 = recMu2.pt();
1323             double eta2 = recMu2.eta();
1324             // This is to avoid nan
1325             if (diffMass == diffMass) {
1326               // Mass relative difference vs Pt and Eta. To be used to extract the true mass resolution
1327               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
1328               mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
1329               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
1330               mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
1331               // This is used for the covariance comparison
1332               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
1333               mapHisto_["hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
1334             } else {
1335               std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
1336             }
1337           }
1338           // Fill with mass resolution from resolution function
1339           double massRes = MuScleFitUtils::massResolution(recMu1, recMu2, MuScleFitUtils::parResol);
1340           mapHisto_["hFunctionResolMass"]->Fill(recMu1, std::pow(massRes, 2), -1);
1341           mapHisto_["hFunctionResolMass"]->Fill(recMu2, std::pow(massRes, 2), +1);
1342         }
1343 
1344         mapHisto_["hMass_P"]->Fill(bestRecRes.mass(), prob);
1345         if (debug_ > 0)
1346           std::cout << "mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
1347         mapHisto_["hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
1348 
1349         mapHisto_["hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
1350         mapHisto_["hMassProbVsMu"]->Fill(recMu1, bestRecRes, -1, prob);
1351         mapHisto_["hMassProbVsMu"]->Fill(recMu2, bestRecRes, +1, prob);
1352         mapHisto_["hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
1353         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu1, bestRecRes, -1, prob);
1354         mapHisto_["hMassProbVsMu_fine"]->Fill(recMu2, bestRecRes, +1, prob);
1355       }
1356     }
1357   }  // end if ResFound
1358 
1359   // Fill the pair
1360   // -------------
1361   if (loopCounter > 0) {
1362     if (debug_ > 0)
1363       std::cout << "[MuScleFit]: filling the pair" << std::endl;
1364     MuScleFitUtils::SavedPair[iev] = std::make_pair(recMu1, recMu2);
1365   }
1366 
1367   iev++;
1368   MuScleFitUtils::iev_++;
1369 
1370   // return kContinue;
1371 }
1372 
1373 bool MuScleFit::checkDeltaR(reco::Particle::LorentzVector& genMu, reco::Particle::LorentzVector& recMu) {
1374   //first is always mu-, second is always mu+
1375   double deltaR =
1376       sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(), genMu.Phi()) +
1377            ((recMu.Eta() - genMu.Eta()) * (recMu.Eta() - genMu.Eta())));
1378   if (deltaR < 0.01)
1379     return true;
1380   else if (debug_ > 0) {
1381     std::cout << "Reco muon " << recMu << " with eta " << recMu.Eta() << " and phi " << recMu.Phi() << std::endl
1382               << " DOES NOT MATCH with generated muon from resonance: " << std::endl
1383               << genMu << " with eta " << genMu.Eta() << " and phi " << genMu.Phi() << std::endl;
1384   }
1385   return false;
1386 }
1387 
1388 void MuScleFit::fillComparisonHistograms(const reco::Particle::LorentzVector& genMu,
1389                                          const reco::Particle::LorentzVector& recMu,
1390                                          const std::string& inputName,
1391                                          const int charge) {
1392   std::string name(inputName + "VSMu");
1393   mapHisto_["hResolPt" + name]->Fill(recMu, (-genMu.Pt() + recMu.Pt()) / genMu.Pt(), charge);
1394   mapHisto_["hResolTheta" + name]->Fill(recMu, (-genMu.Theta() + recMu.Theta()), charge);
1395   mapHisto_["hResolCotgTheta" + name]->Fill(
1396       recMu, (-cos(genMu.Theta()) / sin(genMu.Theta()) + cos(recMu.Theta()) / sin(recMu.Theta())), charge);
1397   mapHisto_["hResolEta" + name]->Fill(recMu, (-genMu.Eta() + recMu.Eta()), charge);
1398   mapHisto_["hResolPhi" + name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);
1399 
1400   // Fill only if it was matched to a genMu and this muon is valid
1401   if ((genMu.Pt() != 0) && (recMu.Pt() != 0)) {
1402     mapHisto_["hPtRecoVsPt" + inputName]->Fill(genMu.Pt(), recMu.Pt());
1403   }
1404 }
1405 
1406 void MuScleFit::applySmearing(reco::Particle::LorentzVector& mu) {
1407   if (MuScleFitUtils::SmearType > 0) {
1408     mu = MuScleFitUtils::applySmearing(mu);
1409     if (debug_ > 0)
1410       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after smearing  Pt = " << mu.Pt() << std::endl;
1411   }
1412 }
1413 
1414 void MuScleFit::applyBias(reco::Particle::LorentzVector& mu, const int charge) {
1415   if (MuScleFitUtils::BiasType > 0) {
1416     mu = MuScleFitUtils::applyBias(mu, charge);
1417     if (debug_ > 0)
1418       std::cout << "Muon #" << MuScleFitUtils::goodmuon << ": after bias      Pt = " << mu.Pt() << std::endl;
1419   }
1420 }
1421 
1422 // Simple method to check parameters consistency. It aborts the job if the parameters
1423 // are not consistent.
1424 void MuScleFit::checkParameters() {
1425   // Fits selection dimension check
1426   if (MuScleFitUtils::doResolFit.size() != maxLoopNumber) {
1427     std::cout << "[MuScleFit-Constructor]: wrong size of resolution fits selector = "
1428               << MuScleFitUtils::doResolFit.size() << std::endl;
1429     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1430     abort();
1431   }
1432   if (MuScleFitUtils::doScaleFit.size() != maxLoopNumber) {
1433     std::cout << "[MuScleFit-Constructor]: wrong size of scale fits selector = " << MuScleFitUtils::doScaleFit.size()
1434               << std::endl;
1435     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1436     abort();
1437   }
1438   if (MuScleFitUtils::doCrossSectionFit.size() != maxLoopNumber) {
1439     std::cout << "[MuScleFit-Constructor]: wrong size of cross section fits selector = "
1440               << MuScleFitUtils::doCrossSectionFit.size() << std::endl;
1441     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1442     abort();
1443   }
1444   if (MuScleFitUtils::doBackgroundFit.size() != maxLoopNumber) {
1445     std::cout << "[MuScleFit-Constructor]: wrong size of background fits selector = "
1446               << MuScleFitUtils::doBackgroundFit.size() << std::endl;
1447     std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
1448     abort();
1449   }
1450 
1451   // Bias parameters: dimension check
1452   // --------------------------------
1453   if ((MuScleFitUtils::BiasType == 1 && MuScleFitUtils::parBias.size() != 2) ||  // linear in pt
1454       (MuScleFitUtils::BiasType == 2 && MuScleFitUtils::parBias.size() != 2) ||  // linear in |eta|
1455       (MuScleFitUtils::BiasType == 3 && MuScleFitUtils::parBias.size() != 4) ||  // sinusoidal in phi
1456       (MuScleFitUtils::BiasType == 4 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and |eta|
1457       (MuScleFitUtils::BiasType == 5 && MuScleFitUtils::parBias.size() != 3) ||  // linear in pt and sinusoidal in phi
1458       (MuScleFitUtils::BiasType == 6 &&
1459        MuScleFitUtils::parBias.size() != 3) ||  // linear in |eta| and sinusoidal in phi
1460       (MuScleFitUtils::BiasType == 7 &&
1461        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and |eta| and sinusoidal in phi
1462       (MuScleFitUtils::BiasType == 8 && MuScleFitUtils::parBias.size() != 4) ||   // linear in pt and parabolic in |eta|
1463       (MuScleFitUtils::BiasType == 9 && MuScleFitUtils::parBias.size() != 2) ||   // exponential in pt
1464       (MuScleFitUtils::BiasType == 10 && MuScleFitUtils::parBias.size() != 3) ||  // parabolic in pt
1465       (MuScleFitUtils::BiasType == 11 &&
1466        MuScleFitUtils::parBias.size() != 4) ||  // linear in pt and sin in phi with chg
1467       (MuScleFitUtils::BiasType == 12 &&
1468        MuScleFitUtils::parBias.size() != 6) ||  // linear in pt and para in plus sin in phi with chg
1469       (MuScleFitUtils::BiasType == 13 &&
1470        MuScleFitUtils::parBias.size() != 8) ||  // linear in pt and para in plus sin in phi with chg
1471       MuScleFitUtils::BiasType < 0 ||
1472       MuScleFitUtils::BiasType > 13) {
1473     std::cout << "[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
1474     abort();
1475   }
1476   // Smear parameters: dimension check
1477   // ---------------------------------
1478   if ((MuScleFitUtils::SmearType == 1 && MuScleFitUtils::parSmear.size() != 3) ||
1479       (MuScleFitUtils::SmearType == 2 && MuScleFitUtils::parSmear.size() != 4) ||
1480       (MuScleFitUtils::SmearType == 3 && MuScleFitUtils::parSmear.size() != 5) ||
1481       (MuScleFitUtils::SmearType == 4 && MuScleFitUtils::parSmear.size() != 6) ||
1482       (MuScleFitUtils::SmearType == 5 && MuScleFitUtils::parSmear.size() != 7) ||
1483       (MuScleFitUtils::SmearType == 6 && MuScleFitUtils::parSmear.size() != 16) ||
1484       (MuScleFitUtils::SmearType == 7 && !MuScleFitUtils::parSmear.empty()) || MuScleFitUtils::SmearType < 0 ||
1485       MuScleFitUtils::SmearType > 7) {
1486     std::cout << "[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
1487     abort();
1488   }
1489   // Protect against bad size of parameters
1490   // --------------------------------------
1491   if (MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolFix.size() ||
1492       MuScleFitUtils::parResol.size() != MuScleFitUtils::parResolOrder.size()) {
1493     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
1494     abort();
1495   }
1496   if (MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleFix.size() ||
1497       MuScleFitUtils::parScale.size() != MuScleFitUtils::parScaleOrder.size()) {
1498     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
1499     abort();
1500   }
1501   if (MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionFix.size() ||
1502       MuScleFitUtils::parCrossSection.size() != MuScleFitUtils::parCrossSectionOrder.size()) {
1503     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1504     abort();
1505   }
1506   if (MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrFix.size() ||
1507       MuScleFitUtils::parBgr.size() != MuScleFitUtils::parBgrOrder.size()) {
1508     std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1509     abort();
1510   }
1511 
1512   // Protect against an incorrect number of resonances
1513   // -------------------------------------------------
1514   if (MuScleFitUtils::resfind.size() != 6) {
1515     std::cout << "[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
1516     abort();
1517   }
1518 }
1519 
1520 bool MuScleFit::selGlobalMuon(const pat::Muon* aMuon) {
1521   reco::TrackRef iTrack = aMuon->innerTrack();
1522   const reco::HitPattern& p = iTrack->hitPattern();
1523 
1524   reco::TrackRef gTrack = aMuon->globalTrack();
1525   const reco::HitPattern& q = gTrack->hitPattern();
1526 
1527   return (  //isMuonInAccept(aMuon) &&// no acceptance cuts!
1528       iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 && q.numberOfValidMuonHits() > 0 &&
1529       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 bool MuScleFit::selTrackerMuon(const pat::Muon* aMuon) {
1536   reco::TrackRef iTrack = aMuon->innerTrack();
1537   const reco::HitPattern& p = iTrack->hitPattern();
1538 
1539   return (  //isMuonInAccept(aMuon) // no acceptance cuts!
1540       iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
1541       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
1542       fabs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
1543       fabs(iTrack->dz()) < 15.0);   //should be done w.r.t. PV!
1544 }
1545 
1546 #include "FWCore/Framework/interface/MakerMacros.h"
1547 DEFINE_FWK_LOOPER(MuScleFit);