Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 /** See header file for a class description
0002  *
0003  *  \author S. Bolognesi - INFN Torino / T. Dorigo, M. De Mattia - INFN Padova
0004  * Revised S. Casasso, E. Migliore - UniTo & INFN Torino
0005  */
0006 // Some notes:
0007 // - M(Z) after simulation needs to be extracted as a function of |y_Z| in order to be
0008 //   a better reference point for calibration. In fact, the variation of PDF with y_Z
0009 //   in production is sizable <---- need to check though.
0010 // - ResHalfWidth needs to be optimized - this depends on the level of background.
0011 // - Background parametrization still to be worked on, so far only a constant (type=1, and
0012 //   parameter 2 fixed to 0) works.
0013 // - weights have to be assigned to dimuon mass values in regions where different resonances
0014 //   overlap, and one has to decide which resonance mass to assign the event to - this until
0015 //   we implement in the fitter a sum of probabilities of an event to belong to different
0016 //   resonances. The weight has to depend on the mass and has relative cross sections of
0017 //   Y(1S), 2S, 3S as parameters. Some overlap is also expected in the J/psi-Psi(2S) region
0018 //   when reconstructing masses with Standalone muons.
0019 //
0020 //   MODS 7/7/08 TD:
0021 //   - changed parametrization of resolution in Pt: from sigma_pt = a*Pt + b*|eta| to
0022 //                                                       sigma_pt = (a*Pt + b*|eta|)*Pt
0023 //                                                  which is more correct (I hope)
0024 //   - changed parametrization of resolution in cotgth: from sigma_cotgth = f(eta) to f(cotgth)
0025 // --------------------------------------------------------------------------------------------
0026 
0027 #include "MuScleFitUtils.h"
0028 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0029 #include "SimDataFormats/Track/interface/SimTrack.h"
0030 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0031 #include "DataFormats/Math/interface/LorentzVector.h"
0032 #include "TString.h"
0033 #include "TFile.h"
0034 #include "TTree.h"
0035 #include "TCanvas.h"
0036 #include "TH2F.h"
0037 #include "TF1.h"
0038 #include "TF2.h"
0039 #include <iostream>
0040 #include <fstream>
0041 #include <memory>  // to use the auto_ptr
0042 
0043 // Includes the definitions of all the bias and scale functions
0044 // These functions are selected in the constructor according
0045 // to the input parameters.
0046 #include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
0047 
0048 // To use callgrind for code profiling uncomment also the following define.
0049 //#define USE_CALLGRIND
0050 #ifdef USE_CALLGRIND
0051 #include "valgrind/callgrind.h"
0052 #endif
0053 
0054 // Lorenzian Peak function
0055 // -----------------------
0056 Double_t lorentzianPeak(Double_t *x, Double_t *par) {
0057   return (0.5 * par[0] * par[1] / TMath::Pi()) /
0058          TMath::Max(1.e-10, (x[0] - par[2]) * (x[0] - par[2]) + .25 * par[1] * par[1]);
0059 }
0060 
0061 // Gaussian function
0062 // -----------------
0063 Double_t Gaussian(Double_t *x, Double_t *par) {
0064   return par[0] * exp(-0.5 * ((x[0] - par[1]) / par[2]) * ((x[0] - par[1]) / par[2]));
0065 }
0066 
0067 // Array with number of parameters in the fitting functions
0068 // (not currently in use)
0069 // --------------------------------------------------------
0070 //const int nparsResol[2] = {6, 4};
0071 //const int nparsScale[13] = {2, 2, 2, 3, 3, 3, 4, 4, 2, 3, 4, 6, 8};
0072 //const int nparsBgr[3] = {1, 2, 3};
0073 
0074 // Quantities used for h-value computation
0075 // ---------------------------------------
0076 double mzsum;
0077 double isum;
0078 double f[11][100];
0079 double g[11][100];
0080 
0081 // Lorentzian convoluted with a gaussian:
0082 // --------------------------------------
0083 TF1 *GL = new TF1(
0084     "GL", "0.5/3.1415926*[0]/(pow(x-[1],2)+pow(0.5*[0],2))*exp(-0.5*pow((x-[2])/[3],2))/([3]*sqrt(6.283185))", 0, 1000);
0085 
0086 TF2 *GL2 = new TF2("GL2",
0087                    "0.5/3.1415926*[0]/(pow(x-[1],2)+pow(0.5*[0],2))*exp(-0.5*pow((x-y)/[2],2))/([2]*sqrt(6.283185))",
0088                    0,
0089                    200,
0090                    0,
0091                    200);
0092 
0093 // // Lorentzian convoluted with a gaussian over a linear background:
0094 // // ---------------------------------------------------------------
0095 // TF1 * GLBL = new TF1 ("GLBL",
0096 //   "0.5/3.1415926*[0]/(pow(x-[1],2)+pow(0.5*[0],2))*exp(-0.5*pow((x-[2])/[3],2))/([3]*sqrt(6.283185))+[4]+[5]*x",
0097 //   0, 1000);
0098 
0099 // // Lorentzian convoluted with a gaussian over an exponential background:
0100 // // ---------------------------------------------------------------
0101 // TF1 * GLBE = new TF1 ("GLBE",
0102 //   "0.5/3.1415926*[0]/(pow(x-[1],2)+pow(0.5*[0],2))*exp(-0.5*pow((x-[2])/[3],2))/([3]*sqrt(6.283185))+exp([4]+[5]*x)",
0103 //   0, 1000);
0104 
0105 std::vector<int> MuScleFitUtils::doResolFit;
0106 std::vector<int> MuScleFitUtils::doScaleFit;
0107 std::vector<int> MuScleFitUtils::doCrossSectionFit;
0108 std::vector<int> MuScleFitUtils::doBackgroundFit;
0109 
0110 int MuScleFitUtils::minuitLoop_ = 0;
0111 TH1D *MuScleFitUtils::likelihoodInLoop_ = nullptr;
0112 TH1D *MuScleFitUtils::signalProb_ = nullptr;
0113 TH1D *MuScleFitUtils::backgroundProb_ = nullptr;
0114 
0115 bool MuScleFitUtils::duringMinos_ = false;
0116 
0117 const int MuScleFitUtils::totalResNum = 6;
0118 
0119 int MuScleFitUtils::SmearType = 0;
0120 smearFunctionBase *MuScleFitUtils::smearFunction = nullptr;
0121 int MuScleFitUtils::BiasType = 0;
0122 // No error, we take functions from the same group for bias and scale.
0123 scaleFunctionBase<std::vector<double> > *MuScleFitUtils::biasFunction = nullptr;
0124 int MuScleFitUtils::ResolFitType = 0;
0125 resolutionFunctionBase<double *> *MuScleFitUtils::resolutionFunction = nullptr;
0126 resolutionFunctionBase<std::vector<double> > *MuScleFitUtils::resolutionFunctionForVec = nullptr;
0127 int MuScleFitUtils::ScaleFitType = 0;
0128 scaleFunctionBase<double *> *MuScleFitUtils::scaleFunction = nullptr;
0129 scaleFunctionBase<std::vector<double> > *MuScleFitUtils::scaleFunctionForVec = nullptr;
0130 int MuScleFitUtils::BgrFitType = 0;
0131 
0132 CrossSectionHandler *MuScleFitUtils::crossSectionHandler;
0133 // const int MuScleFitUtils::backgroundFunctionsRegions = 3;
0134 // backgroundFunctionBase * MuScleFitUtils::backgroundFunctionForRegion[MuScleFitUtils::backgroundFunctionsRegions];
0135 // backgroundFunctionBase * MuScleFitUtils::backgroundFunction[MuScleFitUtils::totalResNum];
0136 BackgroundHandler *MuScleFitUtils::backgroundHandler;
0137 std::vector<double> MuScleFitUtils::parBias;
0138 std::vector<double> MuScleFitUtils::parSmear;
0139 std::vector<double> MuScleFitUtils::parResol;
0140 std::vector<double> MuScleFitUtils::parResolStep;
0141 std::vector<double> MuScleFitUtils::parResolMin;
0142 std::vector<double> MuScleFitUtils::parResolMax;
0143 std::vector<double> MuScleFitUtils::parScale;
0144 std::vector<double> MuScleFitUtils::parScaleStep;
0145 std::vector<double> MuScleFitUtils::parScaleMin;
0146 std::vector<double> MuScleFitUtils::parScaleMax;
0147 std::vector<double> MuScleFitUtils::parCrossSection;
0148 std::vector<double> MuScleFitUtils::parBgr;
0149 std::vector<int> MuScleFitUtils::parResolFix;
0150 std::vector<int> MuScleFitUtils::parScaleFix;
0151 std::vector<int> MuScleFitUtils::parCrossSectionFix;
0152 std::vector<int> MuScleFitUtils::parBgrFix;
0153 std::vector<int> MuScleFitUtils::parResolOrder;
0154 std::vector<int> MuScleFitUtils::parScaleOrder;
0155 std::vector<int> MuScleFitUtils::parCrossSectionOrder;
0156 std::vector<int> MuScleFitUtils::parBgrOrder;
0157 
0158 std::vector<int> MuScleFitUtils::resfind;
0159 int MuScleFitUtils::debug = 0;
0160 
0161 bool MuScleFitUtils::ResFound = false;
0162 int MuScleFitUtils::goodmuon = 0;
0163 int MuScleFitUtils::counter_resprob = 0;
0164 
0165 std::vector<std::vector<double> > MuScleFitUtils::parvalue;
0166 
0167 int MuScleFitUtils::FitStrategy = 1;   // Strategy in likelihood fit (1 or 2)
0168 bool MuScleFitUtils::speedup = false;  // Whether to cut corners (no sim study, fewer histos)
0169 
0170 std::vector<std::pair<lorentzVector, lorentzVector> >
0171     MuScleFitUtils::SavedPair;  // Pairs of reconstructed muons making resonances
0172 std::vector<std::pair<lorentzVector, lorentzVector> >
0173     MuScleFitUtils::ReducedSavedPair;  // Pairs of reconstructed muons making resonances inside smaller windows
0174 std::vector<std::pair<lorentzVector, lorentzVector> >
0175     MuScleFitUtils::genPair;  // Pairs of generated muons making resonances
0176 std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >
0177     MuScleFitUtils::SavedPairMuScleFitMuons;  // Pairs of reconstructed muons making resonances
0178 std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> >
0179     MuScleFitUtils::genMuscleFitPair;  // Pairs of generated muons making resonances
0180 //std::vector<GenMuonPair> MuScleFitUtils::genPairMuScleMuons; // Pairs of generated muons making resonances
0181 std::vector<std::pair<lorentzVector, lorentzVector> >
0182     MuScleFitUtils::simPair;  // Pairs of simulated muons making resonances
0183 
0184 // Smearing parameters
0185 // -------------------
0186 double MuScleFitUtils::x[][10000];
0187 
0188 // Probability matrices and normalization values
0189 // ---------------------------------------------
0190 int MuScleFitUtils::nbins = 1000;
0191 double MuScleFitUtils::GLZValue[][1001][1001];
0192 double MuScleFitUtils::GLZNorm[][1001];
0193 double MuScleFitUtils::GLValue[][1001][1001];
0194 double MuScleFitUtils::GLNorm[][1001];
0195 double MuScleFitUtils::ResMaxSigma[];
0196 
0197 // Masses and widths from PDG 2006, half widths to be revised
0198 // NB in particular, halfwidths have to be made a function of muonType
0199 // -------------------------------------------------------------------
0200 const double MuScleFitUtils::mMu2 = 0.011163612;
0201 const double MuScleFitUtils::muMass = 0.105658;
0202 double MuScleFitUtils::ResHalfWidth[];
0203 double MuScleFitUtils::massWindowHalfWidth[][6];
0204 int MuScleFitUtils::MuonType;
0205 int MuScleFitUtils::MuonTypeForCheckMassWindow;
0206 
0207 double MuScleFitUtils::ResGamma[] = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932};
0208 // ATTENTION:
0209 // This is left because the values are used by the BackgroundHandler to define the center of the regions windows,
0210 // but the values used in the code are read computed using the probability histograms ranges.
0211 // The histograms are read after the initialization of the BackgroundHandler (this can be improved so that
0212 // the background handler too could use the new values).
0213 // At this time the values are consistent.
0214 double MuScleFitUtils::ResMinMass[] = {-99, -99, -99, -99, -99, -99};
0215 double MuScleFitUtils::ResMass[] = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969};
0216 // From Summer08 generator production TWiki: https://twiki.cern.ch/twiki/bin/view/CMS/ProductionSummer2008
0217 // - Z->mumu              1.233 nb
0218 // - Upsilon3S->mumu      0.82  nb
0219 // - Upsilon2S->mumu      6.33  nb
0220 // - Upsilon1S->mumu     13.9   nb
0221 // - Prompt Psi2S->mumu   2.169 nb
0222 // - Prompt J/Psi->mumu 127.2   nb
0223 // double MuScleFitUtils::crossSection[] = {1.233, 0.82, 6.33, 13.9, 2.169, 127.2};
0224 // double MuScleFitUtils::crossSection[] = {1.233, 2.07, 6.33, 13.9, 2.169, 127.2};
0225 
0226 unsigned int MuScleFitUtils::loopCounter = 5;
0227 
0228 // According to the pythia manual, there is only a code for the Upsilon and Upsilon'. It does not distinguish
0229 // between Upsilon(2S) and Upsilon(3S)
0230 const unsigned int MuScleFitUtils::motherPdgIdArray[] = {23, 100553, 100553, 553, 100443, 443};
0231 
0232 // double MuScleFitUtils::leftWindowFactor = 1.;
0233 // double MuScleFitUtils::rightWindowFactor = 1.;
0234 
0235 // double MuScleFitUtils::internalLeftWindowFactor = 1.;
0236 // double MuScleFitUtils::internalRightWindowFactor = 1.;
0237 
0238 // int MuScleFitUtils::backgroundWindowEvents_ = 0;
0239 // int MuScleFitUtils::resonanceWindowEvents_ = 0;
0240 
0241 // double MuScleFitUtils::oldEventsOutInRatio_ = 0.;
0242 
0243 bool MuScleFitUtils::scaleFitNotDone_ = true;
0244 
0245 bool MuScleFitUtils::sherpa_ = false;
0246 
0247 bool MuScleFitUtils::useProbsFile_ = true;
0248 
0249 bool MuScleFitUtils::rapidityBinsForZ_ = true;
0250 
0251 bool MuScleFitUtils::separateRanges_ = true;
0252 double MuScleFitUtils::minMuonPt_ = 0.;
0253 double MuScleFitUtils::maxMuonPt_ = 100000000.;
0254 double MuScleFitUtils::minMuonEtaFirstRange_ = -6.;
0255 double MuScleFitUtils::maxMuonEtaFirstRange_ = 6.;
0256 double MuScleFitUtils::minMuonEtaSecondRange_ = -100.;
0257 double MuScleFitUtils::maxMuonEtaSecondRange_ = 100.;
0258 double MuScleFitUtils::deltaPhiMinCut_ = -100.;
0259 double MuScleFitUtils::deltaPhiMaxCut_ = 100.;
0260 
0261 bool MuScleFitUtils::debugMassResol_;
0262 MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents;
0263 
0264 bool MuScleFitUtils::normalizeLikelihoodByEventNumber_ = true;
0265 TMinuit *MuScleFitUtils::rminPtr_ = nullptr;
0266 double MuScleFitUtils::oldNormalization_ = 0.;
0267 unsigned int MuScleFitUtils::normalizationChanged_ = 0;
0268 
0269 bool MuScleFitUtils::startWithSimplex_;
0270 bool MuScleFitUtils::computeMinosErrors_;
0271 bool MuScleFitUtils::minimumShapePlots_;
0272 
0273 int MuScleFitUtils::iev_ = 0;
0274 ///////////////////////////////////////////////////////////////////////////////////////////////
0275 
0276 // Find the best simulated resonance from a vector of simulated muons (SimTracks)
0277 // and return its decay muons
0278 // ------------------------------------------------------------------------------
0279 std::pair<SimTrack, SimTrack> MuScleFitUtils::findBestSimuRes(const std::vector<SimTrack> &simMuons) {
0280   std::pair<SimTrack, SimTrack> simMuFromBestRes;
0281   double maxprob = -0.1;
0282 
0283   // Double loop on muons
0284   // --------------------
0285   for (std::vector<SimTrack>::const_iterator simMu1 = simMuons.begin(); simMu1 != simMuons.end(); simMu1++) {
0286     for (std::vector<SimTrack>::const_iterator simMu2 = simMu1 + 1; simMu2 != simMuons.end(); simMu2++) {
0287       if (((*simMu1).charge() * (*simMu2).charge()) > 0) {
0288         continue;  // this also gets rid of simMu1==simMu2...
0289       }
0290       // Choose the best resonance using its mass. Check Z, Y(3S,2S,1S), Psi(2S), J/Psi in order
0291       // ---------------------------------------------------------------------------------------
0292       double mcomb = ((*simMu1).momentum() + (*simMu2).momentum()).mass();
0293       double Y = ((*simMu1).momentum() + (*simMu2).momentum()).Rapidity();
0294       for (int ires = 0; ires < 6; ires++) {
0295         if (resfind[ires] > 0) {
0296           double prob = massProb(mcomb, Y, ires, 0.);
0297           if (prob > maxprob) {
0298             simMuFromBestRes.first = (*simMu1);
0299             simMuFromBestRes.second = (*simMu2);
0300             maxprob = prob;
0301           }
0302         }
0303       }
0304     }
0305   }
0306 
0307   // Return most likely combination of muons making a resonance
0308   // ----------------------------------------------------------
0309   return simMuFromBestRes;
0310 }
0311 
0312 // Find the best reconstructed resonance from a collection of reconstructed muons
0313 // (MuonCollection) and return its decay muons
0314 // ------------------------------------------------------------------------------
0315 std::pair<MuScleFitMuon, MuScleFitMuon> MuScleFitUtils::findBestRecoRes(const std::vector<MuScleFitMuon> &muons) {
0316   // NB this routine returns the resonance, but it also sets the ResFound flag, which
0317   // is used in MuScleFit to decide whether to use the event or not.
0318   // --------------------------------------------------------------------------------
0319   if (debug > 0)
0320     std::cout << "In findBestRecoRes" << std::endl;
0321   ResFound = false;
0322   std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes;
0323 
0324   // Choose the best resonance using its mass probability
0325   // ----------------------------------------------------
0326   double maxprob = -0.1;
0327   double minDeltaMass = 999999;
0328   std::pair<MuScleFitMuon, MuScleFitMuon> bestMassMuons;
0329   for (std::vector<MuScleFitMuon>::const_iterator Muon1 = muons.begin(); Muon1 != muons.end(); ++Muon1) {
0330     //rc2010
0331     if (debug > 0)
0332       std::cout << "muon_1_charge:" << (*Muon1).charge() << std::endl;
0333     for (std::vector<MuScleFitMuon>::const_iterator Muon2 = Muon1 + 1; Muon2 != muons.end(); ++Muon2) {
0334       //rc2010
0335       if (debug > 0)
0336         std::cout << "after_2" << std::endl;
0337       if ((((*Muon1).charge()) * ((*Muon2).charge())) > 0) {
0338         continue;  // This also gets rid of Muon1==Muon2...
0339       }
0340       // To allow the selection of ranges at negative and positive eta independently we define two
0341       // ranges of eta: (minMuonEtaFirstRange_, maxMuonEtaFirstRange_) and (minMuonEtaSecondRange_, maxMuonEtaSecondRange_).
0342       double ch1 = (*Muon1).charge();
0343       double ch2 = (*Muon2).charge();
0344       double pt1 = (*Muon1).Pt();
0345       double pt2 = (*Muon2).Pt();
0346       double eta1 = (*Muon1).Eta();
0347       double eta2 = (*Muon2).Eta();
0348       if (pt1 >= minMuonPt_ && pt1 < maxMuonPt_ && pt2 >= minMuonPt_ && pt2 < maxMuonPt_ &&
0349           ((eta1 >= minMuonEtaFirstRange_ && eta1 < maxMuonEtaFirstRange_ && eta2 >= minMuonEtaSecondRange_ &&
0350             eta2 < maxMuonEtaSecondRange_ && ch1 >= ch2  // In the configuration file, MuonOne==MuonPlus
0351             ) ||
0352            (eta1 >= minMuonEtaSecondRange_ && eta1 < maxMuonEtaSecondRange_ && eta2 >= minMuonEtaFirstRange_ &&
0353             eta2 < maxMuonEtaFirstRange_ && ch1 < ch2))) {
0354         double mcomb = ((*Muon1).p4() + (*Muon2).p4()).mass();
0355         double Y = ((*Muon1).p4() + (*Muon2).p4()).Rapidity();
0356         if (debug > 1) {
0357           std::cout << "muon1 " << (*Muon1).p4().Px() << ", " << (*Muon1).p4().Py() << ", " << (*Muon1).p4().Pz()
0358                     << ", " << (*Muon1).p4().E() << ", " << (*Muon1).charge() << std::endl;
0359           std::cout << "muon2 " << (*Muon2).p4().Px() << ", " << (*Muon2).p4().Py() << ", " << (*Muon2).p4().Pz()
0360                     << ", " << (*Muon2).p4().E() << ", " << (*Muon2).charge() << std::endl;
0361           std::cout << "mcomb " << mcomb << std::endl;
0362         }
0363         double massResol = 0.;
0364         if (useProbsFile_) {
0365           massResol = massResolution((*Muon1).p4(), (*Muon2).p4(), parResol);
0366         }
0367         double prob = 0;
0368         for (int ires = 0; ires < 6; ires++) {
0369           if (resfind[ires] > 0) {
0370             if (useProbsFile_) {
0371               prob = massProb(mcomb, Y, ires, massResol);
0372             }
0373             if (prob > maxprob) {
0374               if (ch1 < 0) {  // store first the mu minus and then the mu plus
0375                 recMuFromBestRes.first = (*Muon1);
0376                 recMuFromBestRes.second = (*Muon2);
0377               } else {
0378                 recMuFromBestRes.first = (*Muon2);
0379                 recMuFromBestRes.second = (*Muon1);
0380               }
0381               if (debug > 0)
0382                 std::cout << "muon_1_charge (after swapping):" << recMuFromBestRes.first.charge() << std::endl;
0383               ResFound = true;  // NNBB we accept "resonances" even outside mass bounds
0384               maxprob = prob;
0385             }
0386             // if( ResMass[ires] == 0 ) {
0387             //   std::cout << "Error: ResMass["<<ires<<"] = " << ResMass[ires] << std::endl;
0388             //   exit(1);
0389             // }
0390             double deltaMass = std::abs(mcomb - ResMass[ires]) / ResMass[ires];
0391             if (deltaMass < minDeltaMass) {
0392               bestMassMuons = std::make_pair((*Muon1), (*Muon2));
0393               minDeltaMass = deltaMass;
0394             }
0395           }
0396         }
0397       }
0398     }
0399   }
0400   //If outside mass window (maxprob==0) then take the two muons with best invariant mass
0401   //(anyway they will not be used in the likelihood calculation, only to fill plots)
0402   if (!maxprob) {
0403     if (bestMassMuons.first.charge() < 0) {
0404       recMuFromBestRes.first = bestMassMuons.first;
0405       recMuFromBestRes.second = bestMassMuons.second;
0406     } else {
0407       recMuFromBestRes.second = bestMassMuons.first;
0408       recMuFromBestRes.first = bestMassMuons.second;
0409     }
0410   }
0411   return recMuFromBestRes;
0412 }
0413 
0414 // Resolution smearing function called to worsen muon Pt resolution at start
0415 // -------------------------------------------------------------------------
0416 lorentzVector MuScleFitUtils::applySmearing(const lorentzVector &muon) {
0417   double pt = muon.Pt();
0418   double eta = muon.Eta();
0419   double phi = muon.Phi();
0420   double E = muon.E();
0421 
0422   double y[7] = {};
0423   for (int i = 0; i < SmearType + 3; i++) {
0424     y[i] = x[i][goodmuon % 10000];
0425   }
0426 
0427   // Use the smear function selected in the constructor
0428   smearFunction->smear(pt, eta, phi, y, parSmear);
0429 
0430   if (debug > 9) {
0431     std::cout << "Smearing Pt,eta,phi = " << pt << " " << eta << " " << phi << "; x = ";
0432     for (int i = 0; i < SmearType + 3; i++) {
0433       std::cout << y[i];
0434     }
0435     std::cout << std::endl;
0436   }
0437 
0438   double ptEtaPhiE[4] = {pt, eta, phi, E};
0439   return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
0440 }
0441 
0442 // Biasing function called to modify muon Pt scale at the start.
0443 // -------------------------------------------------------------
0444 lorentzVector MuScleFitUtils::applyBias(const lorentzVector &muon, const int chg) {
0445   double ptEtaPhiE[4] = {muon.Pt(), muon.Eta(), muon.Phi(), muon.E()};
0446 
0447   if (MuScleFitUtils::debug > 1)
0448     std::cout << "pt before bias = " << ptEtaPhiE[0] << std::endl;
0449 
0450   // Use functors (although not with the () operator)
0451   // Note that we always pass pt, eta and phi, but internally only the needed
0452   // values are used.
0453   // The functors used are takend from the same group used for the scaling
0454   // thus the name of the method used is "scale".
0455   ptEtaPhiE[0] = biasFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, MuScleFitUtils::parBias);
0456 
0457   if (MuScleFitUtils::debug > 1)
0458     std::cout << "pt after bias = " << ptEtaPhiE[0] << std::endl;
0459 
0460   return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
0461 }
0462 
0463 // Version of applyScale accepting a std::vector<double> of parameters
0464 // --------------------------------------------------------------
0465 lorentzVector MuScleFitUtils::applyScale(const lorentzVector &muon, const std::vector<double> &parval, const int chg) {
0466   double *p = new double[(int)(parval.size())];
0467   // Replaced by auto_ptr, which handles delete at the end
0468   // std::unique_ptr<double> p(new double[(int)(parval.size())]);
0469   // Removed auto_ptr, check massResolution for an explanation.
0470   int id = 0;
0471   for (std::vector<double>::const_iterator it = parval.begin(); it != parval.end(); ++it, ++id) {
0472     //(&*p)[id] = *it;
0473     // Also ok would be (p.get())[id] = *it;
0474     p[id] = *it;
0475   }
0476   lorentzVector tempScaleVec(applyScale(muon, p, chg));
0477   delete[] p;
0478   return tempScaleVec;
0479 }
0480 
0481 // This is called by the likelihood to "taste" different values for additional corrections
0482 // ---------------------------------------------------------------------------------------
0483 lorentzVector MuScleFitUtils::applyScale(const lorentzVector &muon, double *parval, const int chg) {
0484   double ptEtaPhiE[4] = {muon.Pt(), muon.Eta(), muon.Phi(), muon.E()};
0485   int shift = parResol.size();
0486 
0487   if (MuScleFitUtils::debug > 1)
0488     std::cout << "pt before scale = " << ptEtaPhiE[0] << std::endl;
0489 
0490   // the address of parval[shift] is passed as pointer to double. Internally it is used as a normal array, thus:
0491   // array[0] = parval[shift], array[1] = parval[shift+1], ...
0492   ptEtaPhiE[0] = scaleFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, &(parval[shift]));
0493 
0494   if (MuScleFitUtils::debug > 1)
0495     std::cout << "pt after scale = " << ptEtaPhiE[0] << std::endl;
0496 
0497   return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
0498 }
0499 
0500 // Useful function to convert 4-vector coordinates
0501 // -----------------------------------------------
0502 lorentzVector MuScleFitUtils::fromPtEtaPhiToPxPyPz(const double *ptEtaPhiE) {
0503   double px = ptEtaPhiE[0] * cos(ptEtaPhiE[2]);
0504   double py = ptEtaPhiE[0] * sin(ptEtaPhiE[2]);
0505   double tmp = 2 * atan(exp(-ptEtaPhiE[1]));
0506   double pz = ptEtaPhiE[0] * cos(tmp) / sin(tmp);
0507   double E = sqrt(px * px + py * py + pz * pz + muMass * muMass);
0508 
0509   // lorentzVector corrMu(px,py,pz,E);
0510   // To fix memory leaks, this is to be substituted with
0511   // std::unique_ptr<lorentzVector> corrMu(new lorentzVector(px, py, pz, E));
0512 
0513   return lorentzVector(px, py, pz, E);
0514 }
0515 
0516 // Dimuon mass
0517 // -----------
0518 inline double MuScleFitUtils::invDimuonMass(const lorentzVector &mu1, const lorentzVector &mu2) {
0519   return (mu1 + mu2).mass();
0520 }
0521 
0522 // Mass resolution - version accepting a std::vector<double> parval
0523 // -----------------------------------------------------------
0524 double MuScleFitUtils::massResolution(const lorentzVector &mu1,
0525                                       const lorentzVector &mu2,
0526                                       const std::vector<double> &parval) {
0527   // double * p = new double[(int)(parval.size())];
0528   // Replaced by auto_ptr, which handles delete at the end
0529   // --------- //
0530   // ATTENTION //
0531   // --------- //
0532   // auto_ptr calls delete, not delete[] and thus it must
0533   // not be used with arrays. There are alternatives see
0534   // e.g.: http://www.gotw.ca/gotw/042.htm. The best
0535   // alternative seems to be to switch to vector though.
0536   // std::unique_ptr<double> p(new double[(int)(parval.size())]);
0537 
0538   double *p = new double[(int)(parval.size())];
0539   std::vector<double>::const_iterator it = parval.begin();
0540   int id = 0;
0541   for (; it != parval.end(); ++it, ++id) {
0542     // (&*p)[id] = *it;
0543     p[id] = *it;
0544   }
0545   double massRes = massResolution(mu1, mu2, p);
0546   delete[] p;
0547   return massRes;
0548 }
0549 
0550 /**
0551  * We use the following formula: <br>
0552  *
0553  * M = sqrt ( (E1+E2)^2 - (P1+P2)^2 ) <br>
0554  *
0555  * where we express E and P as a function of Pt, phi, and theta: <br>
0556  *
0557  * E  = sqrt ( Pt^2*(1+cotg(theta)^2) + M_mu^2 ) <br>
0558  * Px = Pt*cos(phi), Py = Pt*sin(phi), Pz = Pt*cotg(theta) <br>
0559  *
0560  * from which we find <br>
0561  *
0562  * M = sqrt( 2*M_mu^2 + 2*sqrt(Pt1^2/sin(theta1)^2 + M_mu^2)*sqrt(Pt2^2/sin(theta2)^2 + M_mu^2) -
0563  *           2*Pt1*Pt2* ( cos(phi1-phi2) + cotg(theta1)*cotg(theta2) ) )
0564  *
0565  * and derive WRT Pt1, Pt2, phi1, phi2, theta1, theta2 to get the resolution.
0566  */
0567 double MuScleFitUtils::massResolution(const lorentzVector &mu1, const lorentzVector &mu2, double *parval) {
0568   double mass = (mu1 + mu2).mass();
0569   double pt1 = mu1.Pt();
0570   double phi1 = mu1.Phi();
0571   double eta1 = mu1.Eta();
0572   double theta1 = 2 * atan(exp(-eta1));
0573   double pt2 = mu2.Pt();
0574   double phi2 = mu2.Phi();
0575   double eta2 = mu2.Eta();
0576   double theta2 = 2 * atan(exp(-eta2));
0577 
0578   double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
0579                        sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0580                    pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
0581                   mass;
0582   double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
0583                        sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0584                    pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
0585                   mass;
0586   double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
0587   double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
0588   double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
0589                            sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0590                        pt1 * pt2 * cos(theta2) / sin(theta2)) /
0591                       mass;
0592   double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
0593                            sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0594                        pt2 * pt1 * cos(theta1) / sin(theta1)) /
0595                       mass;
0596 
0597   if (debugMassResol_) {
0598     massResolComponents.dmdpt1 = dmdpt1;
0599     massResolComponents.dmdpt2 = dmdpt2;
0600     massResolComponents.dmdphi1 = dmdphi1;
0601     massResolComponents.dmdphi2 = dmdphi2;
0602     massResolComponents.dmdcotgth1 = dmdcotgth1;
0603     massResolComponents.dmdcotgth2 = dmdcotgth2;
0604   }
0605 
0606   // Resolution parameters:
0607   // ----------------------
0608   double sigma_pt1 = resolutionFunction->sigmaPt(pt1, eta1, parval);
0609   double sigma_pt2 = resolutionFunction->sigmaPt(pt2, eta2, parval);
0610   double sigma_phi1 = resolutionFunction->sigmaPhi(pt1, eta1, parval);
0611   double sigma_phi2 = resolutionFunction->sigmaPhi(pt2, eta2, parval);
0612   double sigma_cotgth1 = resolutionFunction->sigmaCotgTh(pt1, eta1, parval);
0613   double sigma_cotgth2 = resolutionFunction->sigmaCotgTh(pt2, eta2, parval);
0614   double cov_pt1pt2 = resolutionFunction->covPt1Pt2(pt1, eta1, pt2, eta2, parval);
0615 
0616   // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
0617   // multiply it by pt.
0618   double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
0619                          std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
0620                          std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2) +
0621                          2 * dmdpt1 * dmdpt2 * cov_pt1pt2 * sigma_pt1 * sigma_pt2);
0622 
0623   if (debug > 19) {
0624     std::cout << "  Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1) / sin(theta1) << " - Pt2=" << pt2
0625               << " phi2=" << phi2 << " cotgth2=" << cos(theta2) / sin(theta2) << std::endl;
0626     std::cout << " P[0]=" << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3]
0627               << std::endl;
0628     std::cout << "  Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1=" << sigma_pt1
0629               << " sigma_pt2=" << sigma_pt2 << std::endl;
0630     std::cout << "  Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1=" << sigma_phi1
0631               << " sigma_phi2=" << sigma_phi2 << std::endl;
0632     std::cout << "  Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2 << " sigma_cotgth1=" << sigma_cotgth1
0633               << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
0634     std::cout << "  Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2 << " : " << mass << " +- "
0635               << mass_res << std::endl;
0636   }
0637 
0638   // Debug std::cout
0639   // ----------
0640   bool didit = false;
0641   for (int ires = 0; ires < 6; ires++) {
0642     if (!didit && resfind[ires] > 0 && std::abs(mass - ResMass[ires]) < ResHalfWidth[ires]) {
0643       if (mass_res > ResMaxSigma[ires] && counter_resprob < 100) {
0644         counter_resprob++;
0645         LogDebug("MuScleFitUtils") << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
0646         didit = true;
0647       }
0648     }
0649   }
0650 
0651   return mass_res;
0652 }
0653 
0654 /**
0655  * This method can be used outside MuScleFit. It gets the ResolutionFunction that must have been built with the parameters. <br>
0656  * TO-DO: this method duplicates the code in the previous method. It should be changed to avoid the duplication.
0657  */
0658 double MuScleFitUtils::massResolution(const lorentzVector &mu1,
0659                                       const lorentzVector &mu2,
0660                                       const ResolutionFunction &resolFunc) {
0661   double mass = (mu1 + mu2).mass();
0662   double pt1 = mu1.Pt();
0663   double phi1 = mu1.Phi();
0664   double eta1 = mu1.Eta();
0665   double theta1 = 2 * atan(exp(-eta1));
0666   double pt2 = mu2.Pt();
0667   double phi2 = mu2.Phi();
0668   double eta2 = mu2.Eta();
0669   double theta2 = 2 * atan(exp(-eta2));
0670 
0671   double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
0672                        sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0673                    pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
0674                   mass;
0675   double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
0676                        sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0677                    pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
0678                   mass;
0679   double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
0680   double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
0681   double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
0682                            sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
0683                        pt1 * pt2 * cos(theta2) / sin(theta2)) /
0684                       mass;
0685   double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
0686                            sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
0687                        pt2 * pt1 * cos(theta1) / sin(theta1)) /
0688                       mass;
0689 
0690   // Resolution parameters:
0691   // ----------------------
0692   double sigma_pt1 = resolFunc.sigmaPt(mu1);
0693   double sigma_pt2 = resolFunc.sigmaPt(mu2);
0694   double sigma_phi1 = resolFunc.sigmaPhi(mu1);
0695   double sigma_phi2 = resolFunc.sigmaPhi(mu2);
0696   double sigma_cotgth1 = resolFunc.sigmaCotgTh(mu1);
0697   double sigma_cotgth2 = resolFunc.sigmaCotgTh(mu2);
0698 
0699   // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
0700   // multiply it by pt.
0701   double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
0702                          std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
0703                          std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2));
0704 
0705   return mass_res;
0706 }
0707 
0708 // Mass probability - version with linear background included, accepts std::vector<double> parval
0709 // -----------------------------------------------------------------------------------------
0710 double MuScleFitUtils::massProb(const double &mass,
0711                                 const double &resEta,
0712                                 const double &rapidity,
0713                                 const double &massResol,
0714                                 const std::vector<double> &parval,
0715                                 const bool doUseBkgrWindow,
0716                                 const double &eta1,
0717                                 const double &eta2) {
0718 #ifdef USE_CALLGRIND
0719   CALLGRIND_START_INSTRUMENTATION;
0720 #endif
0721 
0722   double *p = new double[(int)(parval.size())];
0723   // Replaced by auto_ptr, which handles delete at the end
0724   // Removed auto_ptr, check massResolution for an explanation.
0725   // std::unique_ptr<double> p(new double[(int)(parval.size())]);
0726 
0727   std::vector<double>::const_iterator it = parval.begin();
0728   int id = 0;
0729   for (; it != parval.end(); ++it, ++id) {
0730     // (&*p)[id] = *it;
0731     p[id] = *it;
0732   }
0733   // p must be passed by value as below:
0734   double massProbability = massProb(mass, resEta, rapidity, massResol, p, doUseBkgrWindow, eta1, eta2);
0735   delete[] p;
0736 
0737 #ifdef USE_CALLGRIND
0738   CALLGRIND_STOP_INSTRUMENTATION;
0739   CALLGRIND_DUMP_STATS;
0740 #endif
0741 
0742   return massProbability;
0743 }
0744 
0745 /**
0746  * After the introduction of the rapidity bins for the Z the probability method works in the following way:
0747  * - if passing iRes == 0, iY is used to select the rapidity bin
0748  * - if passing iRes != 0, iY is used to select the resonance
0749  */
0750 double MuScleFitUtils::probability(const double &mass,
0751                                    const double &massResol,
0752                                    const double GLvalue[][1001][1001],
0753                                    const double GLnorm[][1001],
0754                                    const int iRes,
0755                                    const int iY) {
0756   if (iRes == 0 && iY > 23) {
0757     std::cout << "WARNING: rapidity bin selected = " << iY << " but there are only histograms for the first 24 bins"
0758               << std::endl;
0759   }
0760 
0761   double PS = 0.;
0762   bool insideProbMassWindow = true;
0763   // Interpolate the four values of GLZValue[] in the
0764   // grid square within which the (mass,sigma) values lay
0765   // ----------------------------------------------------
0766   // This must be done with respect to the width used in the computation of the probability distribution,
0767   // so that the bin 0 really matches the bin 0 of that distribution.
0768   //  double fracMass = (mass-(ResMass[iRes]-ResHalfWidth[iRes]))/(2*ResHalfWidth[iRes]);
0769   double fracMass = (mass - ResMinMass[iRes]) / (2 * ResHalfWidth[iRes]);
0770   if (debug > 1)
0771     std::cout << std::setprecision(9) << "mass ResMinMass[iRes] ResHalfWidth[iRes] ResHalfWidth[iRes]" << mass << " "
0772               << ResMinMass[iRes] << " " << ResHalfWidth[iRes] << " " << ResHalfWidth[iRes] << std::endl;
0773   int iMassLeft = (int)(fracMass * (double)nbins);
0774   int iMassRight = iMassLeft + 1;
0775   double fracMassStep = (double)nbins * (fracMass - (double)iMassLeft / (double)nbins);
0776   if (debug > 1)
0777     std::cout << "nbins iMassLeft fracMass " << nbins << " " << iMassLeft << " " << fracMass << std::endl;
0778 
0779   // Simple protections for the time being: the region where we fit should not include
0780   // values outside the boundaries set by ResMass-ResHalfWidth : ResMass+ResHalfWidth
0781   // ---------------------------------------------------------------------------------
0782   if (iMassLeft < 0) {
0783     edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassLeft=" << iMassLeft
0784                                 << "; mass = " << mass << " and bounds are " << ResMinMass[iRes] << ":"
0785                                 << ResMinMass[iRes] + 2 * ResHalfWidth[iRes] << " - iMassLeft set to 0" << std::endl;
0786     iMassLeft = 0;
0787     iMassRight = 1;
0788     insideProbMassWindow = false;
0789   }
0790   if (iMassRight > nbins) {
0791     edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassRight=" << iMassRight
0792                                 << "; mass = " << mass << " and bounds are " << ResMinMass[iRes] << ":"
0793                                 << ResMass[iRes] + 2 * ResHalfWidth[iRes] << " - iMassRight set to " << nbins - 1
0794                                 << std::endl;
0795     iMassLeft = nbins - 1;
0796     iMassRight = nbins;
0797     insideProbMassWindow = false;
0798   }
0799   double fracSigma = (massResol / ResMaxSigma[iRes]);
0800   int iSigmaLeft = (int)(fracSigma * (double)nbins);
0801   int iSigmaRight = iSigmaLeft + 1;
0802   double fracSigmaStep = (double)nbins * (fracSigma - (double)iSigmaLeft / (double)nbins);
0803   // std::cout << "massResol = " << massResol << std::endl;
0804   // std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
0805   // std::cout << "fracSigma = " << fracSigma << std::endl;
0806   // std::cout << "nbins = " << nbins << std::endl;
0807   // std::cout << "ISIGMALEFT = " << iSigmaLeft << std::endl;
0808   // std::cout << "ISIGMARIGHT = " << iSigmaRight << std::endl;
0809   // std::cout << "fracSigmaStep = " << fracSigmaStep << std::endl;
0810 
0811   // Simple protections for the time being: they should not affect convergence, since
0812   // ResMaxSigma is set to very large values, and if massResol exceeds them the fit
0813   // should not get any prize for that (for large sigma, the prob. distr. becomes flat)
0814   // ----------------------------------------------------------------------------------
0815   if (iSigmaLeft < 0) {
0816     edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaLeft=" << iSigmaLeft
0817                                 << ", with massResol = " << massResol
0818                                 << " and ResMaxSigma[iRes] = " << ResMaxSigma[iRes] << " -  iSigmaLeft set to 0"
0819                                 << std::endl;
0820     iSigmaLeft = 0;
0821     iSigmaRight = 1;
0822   }
0823   if (iSigmaRight > nbins) {
0824     if (counter_resprob < 100)
0825       edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaRight=" << iSigmaRight
0826                                   << ", with massResol = " << massResol
0827                                   << " and ResMaxSigma[iRes] = " << ResMaxSigma[iRes] << " -  iSigmaRight set to "
0828                                   << nbins - 1 << std::endl;
0829     iSigmaLeft = nbins - 1;
0830     iSigmaRight = nbins;
0831   }
0832 
0833   // If f11,f12,f21,f22 are the values at the four corners, one finds by linear interpolation the
0834   // formula below for PS
0835   // --------------------------------------------------------------------------------------------
0836   if (insideProbMassWindow) {
0837     double f11 = 0.;
0838     if (GLnorm[iY][iSigmaLeft] > 0)
0839       f11 = GLvalue[iY][iMassLeft][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
0840     double f12 = 0.;
0841     if (GLnorm[iY][iSigmaRight] > 0)
0842       f12 = GLvalue[iY][iMassLeft][iSigmaRight] / GLnorm[iY][iSigmaRight];
0843     double f21 = 0.;
0844     if (GLnorm[iY][iSigmaLeft] > 0)
0845       f21 = GLvalue[iY][iMassRight][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
0846     double f22 = 0.;
0847     if (GLnorm[iY][iSigmaRight] > 0)
0848       f22 = GLvalue[iY][iMassRight][iSigmaRight] / GLnorm[iY][iSigmaRight];
0849     PS = f11 + (f12 - f11) * fracSigmaStep + (f21 - f11) * fracMassStep +
0850          (f22 - f21 - f12 + f11) * fracMassStep * fracSigmaStep;
0851     if (PS > 0.1 || debug > 1)
0852       LogDebug("MuScleFitUtils") << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22=" << f11 << " " << f12 << " "
0853                                  << f21 << " " << f22 << " "
0854                                  << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR=" << iSigmaLeft
0855                                  << " " << iSigmaRight << " GLvalue[" << iY << "][" << iMassLeft
0856                                  << "] = " << GLvalue[iY][iMassLeft][iSigmaLeft] << " GLnorm[" << iY << "]["
0857                                  << iSigmaLeft << "] = " << GLnorm[iY][iSigmaLeft] << std::endl;
0858 
0859     //     if (PS>0.1) std::cout << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
0860     //                      << f11 << " " << f12 << " " << f21 << " " << f22 << " "
0861     //                      << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
0862     //                      << iSigmaLeft << " " << iSigmaRight << " GLV,GLN="
0863     //                      << GLvalue[iY][iMassLeft][iSigmaLeft]
0864     //                      << " " << GLnorm[iY][iSigmaLeft] << std::endl;
0865 
0866   } else {
0867     edm::LogInfo("probability") << "outside mass probability window. Setting PS[" << iRes << "] = 0" << std::endl;
0868   }
0869 
0870   //   if( PS != PS ) {
0871   //     std::cout << "ERROR: PS = " << PS << " for iRes = " << iRes << std::endl;
0872 
0873   //     std::cout << "mass = " << mass << ", massResol = " << massResol << std::endl;
0874   //     std::cout << "fracMass = " << fracMass << ", iMassLeft = " << iMassLeft
0875   //          << ", iMassRight = " << iMassRight << ", fracMassStep = " << fracMassStep << std::endl;
0876   //     std::cout << "fracSigma = " << fracSigma << ", iSigmaLeft = " << iSigmaLeft
0877   //          << ", iSigmaRight = " << iSigmaRight << ", fracSigmaStep = " << fracSigmaStep << std::endl;
0878   //     std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
0879 
0880   //     std::cout << "GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
0881   //          << " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;
0882   //   }
0883 
0884   return PS;
0885 }
0886 
0887 // Mass probability - version with linear background included
0888 // ----------------------------------------------------------
0889 double MuScleFitUtils::massProb(const double &mass,
0890                                 const double &resEta,
0891                                 const double &rapidity,
0892                                 const double &massResol,
0893                                 double *parval,
0894                                 const bool doUseBkgrWindow,
0895                                 const double &eta1,
0896                                 const double &eta2) {
0897   // This routine computes the likelihood that a given measured mass "measMass" is
0898   // the result of a reference mass ResMass[] if the resolution
0899   // expected for the two muons is massResol.
0900   // This version includes two parameters (the last two in parval, by default)
0901   // to size up the background fraction and its relative normalization with respect
0902   // to the signal shape.
0903   //
0904   // We model the signal probability with a Lorentz L(M,H) of resonance mass M and natural width H
0905   // convoluted with a gaussian G(m,s) of measured mass m and expected mass resolution s,
0906   // by integrating over the intersection of the supports of L and G (which can be made to coincide with
0907   // the region where L is non-zero, given that H<<s in most cases) the product L(M,H)*G(m,s) dx as follows:
0908   //
0909   //   GL(m,s) = Int(M-10H,M+10H) [ L(x-M,H) * G(x-m,s) ] dx
0910   //
0911   // The above convolution is computed numerically by an independent root macro, Probs.C, which outputs
0912   // the values in six 1001x1001 grids, one per resonance.
0913   //
0914   // NB THe following block of explanations for background models is outdated, see detailed
0915   // explanations where the code computes PB.
0916   // +++++++++++++++++++++++
0917   //   For the background, instead, we have two choices: a linear and an exponential model.
0918   //     * For the linear model, we choose a one-parameter form whereby the line is automatically normalized
0919   //       in the support [x1,x2] where we defined our "signal region", as follows:
0920   //
0921   //         B(x;b) = 1/(x2-x1) + {x - (x2+x1)/2} * b
0922   //
0923   //       Defined as above, B(x) is a line passing for the point of coordinates (x1+x2)/2, 1/(x2-x1),
0924   //       whose slope b has as support the interval ( -2/[(x1-x2)*(x1+x2)], 2/[(x1-x2)*(x1+x2)] )
0925   //       so that B(x) is always positive-definite in [x1,x2].
0926   //
0927   //     * For the exponential model, we define B(x;b) as
0928   //
0929   //         B(x;b) = b * { exp(-b*x) / [exp(-b*x1)-exp(-b*x2)] }
0930   //
0931   //       This one-parameter definition is automatically normalized to unity in [x1,x2], with a parameter
0932   //       b which has to be positive in order for the slope to be negative.
0933   //       Please note that this model is not useful in most circumstances; a more useful form would be one
0934   //       which included a linear component.
0935   // ++++++++++++++++++++++
0936   //
0937   // Once GL(m,s) and B(x;b) are defined, we introduce a further parameter a, such that we can have the
0938   // likelihood control the relative fraction of signal and background. We first normalize GL(m,s) for
0939   // any given s by taking the integral
0940   //
0941   //   Int(x1,x2) GL(m,s) dm = K_s
0942   //
0943   // We then define the probability as
0944   //
0945   //   P(m,s,a,b) = GL(m,s)/K_s * a  +  B(x,b) * (1-a)
0946   //
0947   // with a taking on values in the interval [0,1].
0948   // Defined as above, the probability is well-behaved, in the sense that it has a value between 0 and 1,
0949   // and the four parameters m,s,a,b fully control its shape.
0950   //
0951   // It is to be noted that the formulation above requires the computation of two rather time-consuming
0952   // integrals. The one defining GL(m,s) can be stored in a TH2D and loaded by the constructor from a
0953   // file suitably prepared, and this will save loads of computing time.
0954   // ----------------------------------------------------------------------------------------------------
0955 
0956   double P = 0.;
0957   int crossSectionParShift = parResol.size() + parScale.size();
0958   // Take the relative cross sections
0959   std::vector<double> relativeCrossSections =
0960       crossSectionHandler->relativeCrossSections(&(parval[crossSectionParShift]), resfind);
0961   //   for( unsigned int i=0; i<relativeCrossSections.size(); ++i ) {
0962   //     std::cout << "relativeCrossSections["<<i<<"] = " << relativeCrossSections[i] << std::endl;
0963   //     std::cout << "parval["<<crossSectionParShift+i<<"] = " << parval[crossSectionParShift+i] << std::endl;
0964   //   }
0965 
0966   // int bgrParShift = crossSectionParShift + parCrossSection.size();
0967   int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
0968   double Bgrp1 = 0.;
0969   //   double Bgrp2 = 0.;
0970   //   double Bgrp3 = 0.;
0971 
0972   // NB defined as below, P is a non-rigorous "probability" that we observe
0973   // a dimuon mass "mass", given ResMass[], gamma, and massResol. It is what we need for the
0974   // fit which finds the best resolution parameters, though. A definition which is
0975   // more properly a probability is given below (in massProb2()), where the result
0976   // cannot be used to fit resolution parameters because the fit would always prefer
0977   // to set the res parameters to the minimum possible value (best resolution),
0978   // to have a probability close to one of observing any mass.
0979   // -------------------------------------------------------------------------------
0980 
0981   // Determine what resonance(s) we have to deal with
0982   // NB for now we assume equal xs for each resonance
0983   // so we do not assign them different weights
0984   // ------------------------------------------------
0985   double PS[6] = {0.};
0986   double PB = 0.;
0987   double PStot[6] = {0.};
0988 
0989   // Should be removed because it is not used
0990   bool resConsidered[6] = {false};
0991 
0992   bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
0993   // bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
0994 
0995   // First check the Z, which is divided in 24 rapidity bins
0996   // NB max value of Z rapidity to be considered is 2.4 here
0997   // -------------------------------------------------------
0998 
0999   // Do this only if we want to use the rapidity bins for the Z
1000   if (MuScleFitUtils::rapidityBinsForZ_) {
1001     // ATTENTION: cut on Z rapidity at 2.4 since we only have histograms up to that value
1002     // std::pair<double, double> windowFactors = backgroundHandler->windowFactors( useBackgroundWindow, 0 );
1003     std::pair<double, double> windowBorders = backgroundHandler->windowBorders(useBackgroundWindow, 0);
1004     if (resfind[0] > 0
1005         // && checkMassWindow( mass, 0,
1006         //                     backgroundHandler->resMass( useBackgroundWindow, 0 ),
1007         //                     windowFactors.first, windowFactors.second )
1008         && checkMassWindow(mass, windowBorders.first, windowBorders.second)
1009         // && std::abs(rapidity)<2.4
1010     ) {
1011       int iY = (int)(std::abs(rapidity) * 10.);
1012       if (iY > 23)
1013         iY = 23;
1014 
1015       if (MuScleFitUtils::debug > 1)
1016         std::cout << "massProb:resFound = 0, rapidity bin =" << iY << std::endl;
1017 
1018       // In this case the last value is the rapidity bin
1019       PS[0] = probability(mass, massResol, GLZValue, GLZNorm, 0, iY);
1020 
1021       if (PS[0] != PS[0]) {
1022         std::cout << "ERROR: PS[0] = nan, setting it to 0" << std::endl;
1023         PS[0] = 0;
1024       }
1025 
1026       // std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
1027       //                                           &(parval[bgrParShift]), MuScleFitUtils::totalResNum, 0,
1028       //                                           resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
1029 
1030       std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction(doBackgroundFit[loopCounter],
1031                                                                                   &(parval[bgrParShift]),
1032                                                                                   MuScleFitUtils::totalResNum,
1033                                                                                   0,
1034                                                                                   resConsidered,
1035                                                                                   ResMass,
1036                                                                                   ResHalfWidth,
1037                                                                                   MuonType,
1038                                                                                   mass,
1039                                                                                   eta1,
1040                                                                                   eta2);
1041 
1042       Bgrp1 = bgrResult.first;
1043       // When fitting the background we have only one Bgrp1
1044       // When not fitting the background we have many only in a superposition region and this case is treated
1045       // separately after this loop
1046       PB = bgrResult.second;
1047       if (PB != PB)
1048         PB = 0;
1049       PStot[0] = (1 - Bgrp1) * PS[0] + Bgrp1 * PB;
1050 
1051       // PStot[0] *= crossSectionHandler->crossSection(0);
1052       // PStot[0] *= parval[crossSectionParShift];
1053       PStot[0] *= relativeCrossSections[0];
1054       if (MuScleFitUtils::debug > 0)
1055         std::cout << "PStot[" << 0 << "] = "
1056                   << "(1-" << Bgrp1 << ")*" << PS[0] << " + " << Bgrp1 << "*" << PB << " = " << PStot[0]
1057                   << " (reletiveXS = )" << relativeCrossSections[0] << std::endl;
1058     } else {
1059       if (debug > 0) {
1060         std::cout << "Mass = " << mass << " outside range with rapidity = " << rapidity << std::endl;
1061         std::cout << "with resMass = " << backgroundHandler->resMass(useBackgroundWindow, 0)
1062                   << " and left border = " << windowBorders.first << " right border = " << windowBorders.second
1063                   << std::endl;
1064       }
1065     }
1066   }
1067   // Next check the other resonances
1068   // -------------------------------
1069   int firstRes = 1;
1070   if (!MuScleFitUtils::rapidityBinsForZ_)
1071     firstRes = 0;
1072   for (int ires = firstRes; ires < 6; ++ires) {
1073     if (resfind[ires] > 0) {
1074       // First is left, second is right (returns (1,1) in the case of resonances, it could be improved avoiding the call in this case)
1075       // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
1076       std::pair<double, double> windowBorder = backgroundHandler->windowBorders(useBackgroundWindow, ires);
1077       if (checkMassWindow(mass, windowBorder.first, windowBorder.second)) {
1078         if (MuScleFitUtils::debug > 1)
1079           std::cout << "massProb:resFound = " << ires << std::endl;
1080 
1081         // In this case the rapidity value is instead the resonance index again.
1082         PS[ires] = probability(mass, massResol, GLValue, GLNorm, ires, ires);
1083 
1084         std::pair<double, double> bgrResult =
1085             backgroundHandler->backgroundFunction(doBackgroundFit[loopCounter],
1086                                                   &(parval[bgrParShift]),
1087                                                   MuScleFitUtils::totalResNum,
1088                                                   ires,
1089                                                   // resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
1090                                                   resConsidered,
1091                                                   ResMass,
1092                                                   ResHalfWidth,
1093                                                   MuonType,
1094                                                   mass,
1095                                                   eta1,
1096                                                   eta2);
1097         Bgrp1 = bgrResult.first;
1098         PB = bgrResult.second;
1099 
1100         if (PB != PB)
1101           PB = 0;
1102         PStot[ires] = (1 - Bgrp1) * PS[ires] + Bgrp1 * PB;
1103         if (MuScleFitUtils::debug > 0)
1104           std::cout << "PStot[" << ires << "] = "
1105                     << "(1-" << Bgrp1 << ")*" << PS[ires] << " + " << Bgrp1 << "*" << PB << " = " << PStot[ires]
1106                     << std::endl;
1107 
1108         PStot[ires] *= relativeCrossSections[ires];
1109       }
1110     }
1111   }
1112 
1113   for (int i = 0; i < 6; ++i) {
1114     P += PStot[i];
1115   }
1116 
1117   if (MuScleFitUtils::signalProb_ != nullptr && MuScleFitUtils::backgroundProb_ != nullptr) {
1118     double PStotTemp = 0.;
1119     for (int i = 0; i < 6; ++i) {
1120       PStotTemp += PS[i] * relativeCrossSections[i];
1121     }
1122     if (PStotTemp != PStotTemp) {
1123       std::cout << "ERROR: PStotTemp = nan!!!!!!!!!" << std::endl;
1124       int parnumber = (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size());
1125       for (int i = 0; i < 6; ++i) {
1126         std::cout << "PS[i] = " << PS[i] << std::endl;
1127         if (PS[i] != PS[i]) {
1128           std::cout << "mass = " << mass << std::endl;
1129           std::cout << "massResol = " << massResol << std::endl;
1130           for (int j = 0; j < parnumber; ++j) {
1131             std::cout << "parval[" << j << "] = " << parval[j] << std::endl;
1132           }
1133         }
1134       }
1135     }
1136     if (PStotTemp == PStotTemp) {
1137       MuScleFitUtils::signalProb_->SetBinContent(
1138           MuScleFitUtils::minuitLoop_,
1139           MuScleFitUtils::signalProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PStotTemp);
1140     }
1141     if (debug > 0)
1142       std::cout << "mass = " << mass << ", P = " << P << ", PStot = " << PStotTemp << ", PB = " << PB
1143                 << ", bgrp1 = " << Bgrp1 << std::endl;
1144 
1145     MuScleFitUtils::backgroundProb_->SetBinContent(
1146         MuScleFitUtils::minuitLoop_, MuScleFitUtils::backgroundProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PB);
1147   }
1148   return P;
1149 }
1150 
1151 // Method to check if the mass value is within the mass window of the i-th resonance.
1152 // inline bool MuScleFitUtils::checkMassWindow( const double & mass, const int ires, const double & resMass, const double & leftFactor, const double & rightFactor )
1153 // {
1154 //   return( mass-resMass > -leftFactor*massWindowHalfWidth[MuonTypeForCheckMassWindow][ires]
1155 //           && mass-resMass < rightFactor*massWindowHalfWidth[MuonTypeForCheckMassWindow][ires] );
1156 // }
1157 inline bool MuScleFitUtils::checkMassWindow(const double &mass, const double &leftBorder, const double &rightBorder) {
1158   return ((mass > leftBorder) && (mass < rightBorder));
1159 }
1160 
1161 // Function that returns the weight for a muon pair
1162 // ------------------------------------------------
1163 double MuScleFitUtils::computeWeight(const double &mass, const int iev, const bool doUseBkgrWindow) {
1164   // Compute weight for this event
1165   // -----------------------------
1166   double weight = 0.;
1167 
1168   // Take the highest-mass resonance within bounds
1169   // NB this must be revised once credible estimates of the relative xs of Y(1S), (2S), and (3S)
1170   // are made. Those are priors in the decision of which resonance to assign to an in-between event.
1171   // -----------------------------------------------------------------------------------------------
1172 
1173   if (doUseBkgrWindow && (debug > 0))
1174     std::cout << "using backgrond window for mass = " << mass << std::endl;
1175   // bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
1176   bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
1177 
1178   for (int ires = 0; ires < 6; ires++) {
1179     if (resfind[ires] > 0 && weight == 0.) {
1180       // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
1181       std::pair<double, double> windowBorder = backgroundHandler->windowBorders(useBackgroundWindow, ires);
1182       // if( checkMassWindow(mass, ires, backgroundHandler->resMass( useBackgroundWindow, ires ),
1183       //                     windowFactor.first, windowFactor.second) ) {
1184       if (checkMassWindow(mass, windowBorder.first, windowBorder.second)) {
1185         weight = 1.0;
1186         if (doUseBkgrWindow && (debug > 0))
1187           std::cout << "setting weight to = " << weight << std::endl;
1188       }
1189     }
1190   }
1191 
1192   return weight;
1193 }
1194 
1195 // Likelihood minimization routine
1196 // -------------------------------
1197 void MuScleFitUtils::minimizeLikelihood() {
1198   // Output file with fit parameters resulting from minimization
1199   // -----------------------------------------------------------
1200   std::ofstream FitParametersFile;
1201   FitParametersFile.open("FitParameters.txt", std::ios::app);
1202   FitParametersFile << "Fitting with resolution, scale, bgr function # " << ResolFitType << " " << ScaleFitType << " "
1203                     << BgrFitType << " - Iteration " << loopCounter << std::endl;
1204 
1205   // Fill parvalue and other vectors needed for the fitting
1206   // ------------------------------------------------------
1207 
1208   // ----- //
1209   // FIXME //
1210   // ----- //
1211   // this was changed to verify the possibility that fixed parameters influence the errors.
1212   // It must be 0 otherwise the parameters for resonances will not be passed by minuit (will be always 0).
1213   // Should be removed.
1214   int parForResonanceWindows = 0;
1215   // int parnumber = (int)(parResol.size()+parScale.size()+parCrossSection.size()+parBgr.size() - parForResonanceWindows);
1216   int parnumber =
1217       (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size() - parForResonanceWindows);
1218 
1219   int parnumberAll = (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size());
1220 
1221   // parvalue is a std::vector<std::vector<double> > storing all the parameters from all the loops
1222   parvalue.push_back(parResol);
1223   std::vector<double> *tmpVec = &(parvalue.back());
1224 
1225   // If this is not the first loop we want to start from neutral values
1226   // Otherwise the scale will start with values correcting again a bias
1227   // that is already corrected.
1228   if (scaleFitNotDone_) {
1229     tmpVec->insert(tmpVec->end(), parScale.begin(), parScale.end());
1230     std::cout << "scaleFitNotDone: tmpVec->size() = " << tmpVec->size() << std::endl;
1231   } else {
1232     scaleFunction->resetParameters(tmpVec);
1233     std::cout << "scaleFitDone: tmpVec->size() = " << tmpVec->size() << std::endl;
1234   }
1235   tmpVec->insert(tmpVec->end(), parCrossSection.begin(), parCrossSection.end());
1236   tmpVec->insert(tmpVec->end(), parBgr.begin(), parBgr.end());
1237   int i = 0;
1238   std::vector<double>::const_iterator it = tmpVec->begin();
1239   for (; it != tmpVec->end(); ++it, ++i) {
1240     std::cout << "tmpVec[" << i << "] = " << *it << std::endl;
1241   }
1242 
1243   // Empty vector of size = number of cross section fitted parameters. Note that the cross section
1244   // fit works in a different way than the others and it uses ratios of the paramters passed via cfg.
1245   // We use this empty vector for compatibility with the rest of the structure.
1246   std::vector<int> crossSectionParNumSizeVec(MuScleFitUtils::crossSectionHandler->parNum(), 0);
1247 
1248   std::vector<int> parfix(parResolFix);
1249   parfix.insert(parfix.end(), parScaleFix.begin(), parScaleFix.end());
1250   parfix.insert(parfix.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end());
1251   parfix.insert(parfix.end(), parBgrFix.begin(), parBgrFix.end());
1252 
1253   std::vector<int> parorder(parResolOrder);
1254   parorder.insert(parorder.end(), parScaleOrder.begin(), parScaleOrder.end());
1255   parorder.insert(parorder.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end());
1256   parorder.insert(parorder.end(), parBgrOrder.begin(), parBgrOrder.end());
1257 
1258   // This is filled later
1259   std::vector<double> parerr(3 * parnumberAll, 0.);
1260 
1261   if (debug > 19) {
1262     std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters before likelihood " << std::endl;
1263     for (unsigned int i = 0; i < (unsigned int)parnumberAll; i++) {
1264       std::cout << "  Par # " << i << " = " << parvalue[loopCounter][i] << " : free = " << parfix[i]
1265                 << "; order = " << parorder[i] << std::endl;
1266     }
1267   }
1268 
1269   //   // Background rescaling from regions to resonances
1270   //   // -----------------------------------------------
1271   //   // If we are in a loop > 0 and we are not fitting the background, but we have fitted it in the previous iteration
1272   //   if( loopCounter > 0 && !(doBackgroundFit[loopCounter]) && doBackgroundFit[loopCounter-1] ) {
1273   //     // This rescales from regions to resonances
1274   //     int localMuonType = MuonType;
1275   //     if( MuonType > 2 ) localMuonType = 2;
1276   //     backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
1277   //                                 MuScleFitUtils::SavedPair);
1278   //   }
1279 
1280   // Init Minuit
1281   // -----------
1282   TMinuit rmin(parnumber);
1283   rminPtr_ = &rmin;
1284   rmin.SetFCN(likelihood);  // Unbinned likelihood
1285   // Standard initialization of minuit parameters:
1286   // sets input to be $stdin, output to be $stdout
1287   // and saving to a file.
1288   rmin.mninit(5, 6, 7);
1289   int ierror = 0;
1290   int istat;
1291   double arglis[4];
1292   arglis[0] = FitStrategy;  // Strategy 1 or 2
1293   // 1 standard
1294   // 2 try to improve minimum (slower)
1295   rmin.mnexcm("SET STR", arglis, 1, ierror);
1296 
1297   arglis[0] = 10001;
1298   // Set the random seed for the generator used in SEEk to a fixed value for reproducibility
1299   rmin.mnexcm("SET RAN", arglis, 1, ierror);
1300 
1301   // Set fit parameters
1302   // ------------------
1303   double *Start = new double[parnumberAll];
1304   double *Step = new double[parnumberAll];
1305   double *Mini = new double[parnumberAll];
1306   double *Maxi = new double[parnumberAll];
1307   int *ind = new int[parnumberAll];  // Order of release of parameters
1308   TString *parname = new TString[parnumberAll];
1309 
1310   if (!parResolStep.empty() && !parResolMin.empty() && !parResolMax.empty()) {
1311     MuScleFitUtils::resolutionFunctionForVec->setParameters(
1312         Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parResolStep, parResolMin, parResolMax, MuonType);
1313   } else {
1314     MuScleFitUtils::resolutionFunctionForVec->setParameters(
1315         Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, MuonType);
1316   }
1317 
1318   // Take the number of parameters in the resolutionFunction and displace the arrays passed to the scaleFunction
1319   int resParNum = MuScleFitUtils::resolutionFunctionForVec->parNum();
1320 
1321   if (!parScaleStep.empty() && !parScaleMin.empty() && !parScaleMax.empty()) {
1322     MuScleFitUtils::scaleFunctionForVec->setParameters(&(Start[resParNum]),
1323                                                        &(Step[resParNum]),
1324                                                        &(Mini[resParNum]),
1325                                                        &(Maxi[resParNum]),
1326                                                        &(ind[resParNum]),
1327                                                        &(parname[resParNum]),
1328                                                        parScale,
1329                                                        parScaleOrder,
1330                                                        parScaleStep,
1331                                                        parScaleMin,
1332                                                        parScaleMax,
1333                                                        MuonType);
1334   } else {
1335     MuScleFitUtils::scaleFunctionForVec->setParameters(&(Start[resParNum]),
1336                                                        &(Step[resParNum]),
1337                                                        &(Mini[resParNum]),
1338                                                        &(Maxi[resParNum]),
1339                                                        &(ind[resParNum]),
1340                                                        &(parname[resParNum]),
1341                                                        parScale,
1342                                                        parScaleOrder,
1343                                                        MuonType);
1344   }
1345 
1346   // Initialize cross section parameters
1347   int crossSectionParShift = resParNum + MuScleFitUtils::scaleFunctionForVec->parNum();
1348   MuScleFitUtils::crossSectionHandler->setParameters(&(Start[crossSectionParShift]),
1349                                                      &(Step[crossSectionParShift]),
1350                                                      &(Mini[crossSectionParShift]),
1351                                                      &(Maxi[crossSectionParShift]),
1352                                                      &(ind[crossSectionParShift]),
1353                                                      &(parname[crossSectionParShift]),
1354                                                      parCrossSection,
1355                                                      parCrossSectionOrder,
1356                                                      resfind);
1357 
1358   // Initialize background parameters
1359   int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
1360   MuScleFitUtils::backgroundHandler->setParameters(&(Start[bgrParShift]),
1361                                                    &(Step[bgrParShift]),
1362                                                    &(Mini[bgrParShift]),
1363                                                    &(Maxi[bgrParShift]),
1364                                                    &(ind[bgrParShift]),
1365                                                    &(parname[bgrParShift]),
1366                                                    parBgr,
1367                                                    parBgrOrder,
1368                                                    MuonType);
1369 
1370   for (int ipar = 0; ipar < parnumber; ++ipar) {
1371     std::cout << "parname[" << ipar << "] = " << parname[ipar] << std::endl;
1372     std::cout << "Start[" << ipar << "] = " << Start[ipar] << std::endl;
1373     std::cout << "Step[" << ipar << "] = " << Step[ipar] << std::endl;
1374     std::cout << "Mini[" << ipar << "] = " << Mini[ipar] << std::endl;
1375     std::cout << "Maxi[" << ipar << "] = " << Maxi[ipar] << std::endl;
1376 
1377     rmin.mnparm(ipar, parname[ipar], Start[ipar], Step[ipar], Mini[ipar], Maxi[ipar], ierror);
1378 
1379     // Testing without limits
1380     // rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], 0, 0, ierror );
1381   }
1382 
1383   // Do minimization
1384   // ---------------
1385   if (debug > 19)
1386     std::cout << "[MuScleFitUtils-minimizeLikelihood]: Starting minimization" << std::endl;
1387   double fmin;
1388   double fdem;
1389   double errdef;
1390   int npari;
1391   int nparx;
1392   rmin.mnexcm("CALL FCN", arglis, 1, ierror);
1393 
1394   // First, fix all parameters
1395   // -------------------------
1396   if (debug > 19)
1397     std::cout << "[MuScleFitUtils-minimizeLikelihood]: First fix all parameters ...";
1398   for (int ipar = 0; ipar < parnumber; ipar++) {
1399     rmin.FixParameter(ipar);
1400   }
1401 
1402   // Then release them in the specified order and refit
1403   // --------------------------------------------------
1404   if (debug > 19)
1405     std::cout << " Then release them in order..." << std::endl;
1406 
1407   TString name;
1408   double pval;
1409   double pmin;
1410   double pmax;
1411   double errp;
1412   double errl;
1413   double errh;
1414   int ivar;
1415   double erro;
1416   double cglo;
1417   int n_times = 0;
1418   // n_times = number of loops required to unlock all parameters.
1419 
1420   if (debug > 19)
1421     std::cout << "Before scale parNum" << std::endl;
1422   int scaleParNum = scaleFunction->parNum();
1423   if (debug > 19)
1424     std::cout << "After scale parNum" << std::endl;
1425   // edm::LogInfo("minimizeLikelihood") << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
1426   // edm::LogInfo("minimizeLikelihood") << "number of parameters for resolutionFunction = " << resParNum << std::endl;
1427   // edm::LogInfo("minimizeLikelihood") << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
1428   // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
1429   std::cout << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
1430   std::cout << "number of parameters for resolutionFunction = " << resParNum << std::endl;
1431   std::cout << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
1432   std::cout << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
1433   // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << backgroundFunction->parNum() << std::endl;
1434 
1435   for (int i = 0; i < parnumber; i++) {
1436     // NB ind[] has been set as parorder[] previously
1437     if (n_times < ind[i]) {
1438       edm::LogInfo("minimizeLikelihood") << "n_times = " << n_times << ", ind[" << i << "] = " << ind[i]
1439                                          << ", scaleParNum = " << scaleParNum << ", doScaleFit[" << loopCounter
1440                                          << "] = " << doScaleFit[loopCounter] << std::endl;
1441       // Set the n_times only if we will do the fit
1442       if (i < resParNum) {
1443         if (doResolFit[loopCounter])
1444           n_times = ind[i];
1445       } else if (i < resParNum + scaleParNum) {
1446         if (doScaleFit[loopCounter])
1447           n_times = ind[i];
1448       } else if (doBackgroundFit[loopCounter])
1449         n_times = ind[i];
1450     }
1451   }
1452   for (int iorder = 0; iorder < n_times + 1; iorder++) {  // Repeat fit n_times times
1453     std::cout << "Starting minimization " << iorder << " of " << n_times << std::endl;
1454 
1455     bool somethingtodo = false;
1456 
1457     // Use parameters from cfg to select which fit to do
1458     // -------------------------------------------------
1459     if (doResolFit[loopCounter]) {
1460       // Release resolution parameters and fit them
1461       // ------------------------------------------
1462       for (unsigned int ipar = 0; ipar < parResol.size(); ++ipar) {
1463         if (parfix[ipar] == 0 && ind[ipar] == iorder) {
1464           rmin.Release(ipar);
1465           somethingtodo = true;
1466         }
1467       }
1468     }
1469     if (doScaleFit[loopCounter]) {
1470       // Release scale parameters and fit them
1471       // -------------------------------------
1472       for (unsigned int ipar = parResol.size(); ipar < parResol.size() + parScale.size(); ++ipar) {
1473         if (parfix[ipar] == 0 && ind[ipar] == iorder) {  // parfix=0 means parameter is free
1474           rmin.Release(ipar);
1475           somethingtodo = true;
1476         }
1477       }
1478       scaleFitNotDone_ = false;
1479     }
1480     unsigned int crossSectionParShift = parResol.size() + parScale.size();
1481     if (doCrossSectionFit[loopCounter]) {
1482       // Release cross section parameters and fit them
1483       // ---------------------------------------------
1484       // Note that only cross sections of resonances that are being fitted are released
1485       bool doCrossSection =
1486           crossSectionHandler->releaseParameters(rmin, resfind, parfix, ind, iorder, crossSectionParShift);
1487       if (doCrossSection)
1488         somethingtodo = true;
1489     }
1490     if (doBackgroundFit[loopCounter]) {
1491       // Release background parameters and fit them
1492       // ------------------------------------------
1493       // for( int ipar=parResol.size()+parScale.size(); ipar<parnumber; ++ipar ) {
1494       // Free only the parameters for the regions, as the resonances intervals are never used to fit the background
1495       unsigned int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
1496       for (unsigned int ipar = bgrParShift; ipar < bgrParShift + backgroundHandler->regionsParNum(); ++ipar) {
1497         // Release only those parameters for the resonances we are fitting
1498         if (parfix[ipar] == 0 && ind[ipar] == iorder &&
1499             backgroundHandler->unlockParameter(resfind, ipar - bgrParShift)) {
1500           rmin.Release(ipar);
1501           somethingtodo = true;
1502         }
1503       }
1504     }
1505 
1506     // OK, now do minimization if some parameter has been released
1507     // -----------------------------------------------------------
1508     if (somethingtodo) {
1509       // #ifdef DEBUG
1510 
1511       std::stringstream fileNum;
1512       fileNum << loopCounter;
1513 
1514       minuitLoop_ = 0;
1515       char name[50];
1516       sprintf(name, "likelihoodInLoop_%d_%d", loopCounter, iorder);
1517       TH1D *tempLikelihoodInLoop = new TH1D(name, "likelihood value in minuit loop", 10000, 0, 10000);
1518       likelihoodInLoop_ = tempLikelihoodInLoop;
1519       char signalProbName[50];
1520       sprintf(signalProbName, "signalProb_%d_%d", loopCounter, iorder);
1521       TH1D *tempSignalProb = new TH1D(signalProbName, "signal probability", 10000, 0, 10000);
1522       signalProb_ = tempSignalProb;
1523       char backgroundProbName[50];
1524       sprintf(backgroundProbName, "backgroundProb_%d_%d", loopCounter, iorder);
1525       TH1D *tempBackgroundProb = new TH1D(backgroundProbName, "background probability", 10000, 0, 10000);
1526       backgroundProb_ = tempBackgroundProb;
1527       // #endif
1528 
1529       // Before we start the minimization we create a list of events with only the events inside a smaller
1530       // window than the one in which the probability is != 0. We will compute the probability for all those
1531       // events and hopefully the margin will avoid them to get a probability = 0 (which causes a discontinuity
1532       // in the likelihood function). The width of this smaller window must be optimized, but we can start using
1533       // an 90% of the normalization window.
1534       double protectionFactor = 0.9;
1535 
1536       MuScleFitUtils::ReducedSavedPair.clear();
1537       for (unsigned int nev = 0; nev < MuScleFitUtils::SavedPair.size(); ++nev) {
1538         const lorentzVector *recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
1539         const lorentzVector *recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
1540         double mass = MuScleFitUtils::invDimuonMass(*recMu1, *recMu2);
1541         // Test all resonances to see if the mass is inside at least one of the windows
1542         bool check = false;
1543         for (int ires = 0; ires < 6; ++ires) {
1544           // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( doBackgroundFit[loopCounter], ires );
1545           std::pair<double, double> windowBorder = backgroundHandler->windowBorders(doBackgroundFit[loopCounter], ires);
1546           // if( resfind[ires] && checkMassWindow( mass, ires, backgroundHandler->resMass( doBackgroundFit[loopCounter], ires ),
1547           //                                       0.9*windowFactor.first, 0.9*windowFactor.second ) ) {
1548           // double resMassValue = backgroundHandler->resMass( doBackgroundFit[loopCounter], ires );
1549           // double windowBorderLeft = resMassValue - protectionFactor*(resMassValue - windowBorder.first);
1550           // double windowBorderRight = resMassValue + protectionFactor*(windowBorder.second - resMassValue);
1551           double windowBorderShift = (windowBorder.second - windowBorder.first) * (1 - protectionFactor) / 2.;
1552           double windowBorderLeft = windowBorder.first + windowBorderShift;
1553           double windowBorderRight = windowBorder.second - windowBorderShift;
1554           if (resfind[ires] && checkMassWindow(mass, windowBorderLeft, windowBorderRight)) {
1555             check = true;
1556           }
1557         }
1558         if (check) {
1559           MuScleFitUtils::ReducedSavedPair.push_back(std::make_pair(*recMu1, *recMu2));
1560         }
1561       }
1562       std::cout << "Fitting with " << MuScleFitUtils::ReducedSavedPair.size() << " events" << std::endl;
1563 
1564       // rmin.SetMaxIterations(500*parnumber);
1565 
1566       //Print some informations
1567       std::cout << "MINUIT is starting the minimization for the iteration number " << loopCounter << std::endl;
1568 
1569       //Try to set iterations
1570       //      rmin.SetMaxIterations(100000);
1571 
1572       std::cout << "maxNumberOfIterations (just set) = " << rmin.GetMaxIterations() << std::endl;
1573 
1574       MuScleFitUtils::normalizationChanged_ = 0;
1575 
1576       // Maximum number of iterations
1577       arglis[0] = 100000;
1578       // tolerance
1579       arglis[1] = 0.1;
1580 
1581       // Run simplex first to get an initial estimate of the minimum
1582       if (startWithSimplex_) {
1583         rmin.mnexcm("SIMPLEX", arglis, 0, ierror);
1584       }
1585 
1586       rmin.mnexcm("MIGRAD", arglis, 2, ierror);
1587 
1588       // #ifdef DEBUG
1589       likelihoodInLoop_->Write();
1590       signalProb_->Write();
1591       backgroundProb_->Write();
1592       delete tempLikelihoodInLoop;
1593       delete tempSignalProb;
1594       delete tempBackgroundProb;
1595       likelihoodInLoop_ = nullptr;
1596       signalProb_ = nullptr;
1597       backgroundProb_ = nullptr;
1598       // #endif
1599 
1600       // Compute again the error matrix
1601       rmin.mnexcm("HESSE", arglis, 0, ierror);
1602 
1603       // Peform minos error analysis.
1604       if (computeMinosErrors_) {
1605         duringMinos_ = true;
1606         rmin.mnexcm("MINOS", arglis, 0, ierror);
1607         duringMinos_ = false;
1608       }
1609 
1610       if (normalizationChanged_ > 1) {
1611         std::cout << "WARNING: normalization changed during fit meaning that events exited from the mass window. This "
1612                      "causes a discontinuity in the likelihood function. Please check the scan of the likelihood as a "
1613                      "function of the parameters to see if there are discontinuities around the minimum."
1614                   << std::endl;
1615       }
1616     }
1617 
1618     // bool notWritten = true;
1619     for (int ipar = 0; ipar < parnumber; ipar++) {
1620       rmin.mnpout(ipar, name, pval, erro, pmin, pmax, ivar);
1621       // Save parameters in parvalue[] vector
1622       // ------------------------------------
1623       if (ierror != 0 && debug > 0) {
1624         std::cout << "[MuScleFitUtils-minimizeLikelihood]: ierror!=0, bogus pars" << std::endl;
1625       }
1626       //     for (int ipar=0; ipar<parnumber; ipar++) {
1627       //       rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
1628       parvalue[loopCounter][ipar] = pval;
1629       //     }
1630 
1631       // int ilax2 = 0;
1632       // Double_t val2pl, val2mi;
1633       // rmin.mnmnot (ipar+1, ilax2, val2pl, val2mi);
1634       rmin.mnerrs(ipar, errh, errl, errp, cglo);
1635 
1636       // Set error on params
1637       // -------------------
1638       if (errp != 0) {
1639         parerr[3 * ipar] = errp;
1640       } else {
1641         parerr[3 * ipar] = (((errh) > (std::abs(errl))) ? (errh) : (std::abs(errl)));
1642       }
1643       parerr[3 * ipar + 1] = errl;
1644       parerr[3 * ipar + 2] = errh;
1645 
1646       if (ipar == 0) {
1647         FitParametersFile << " Resolution fit parameters:" << std::endl;
1648       }
1649       if (ipar == int(parResol.size())) {
1650         FitParametersFile << " Scale fit parameters:" << std::endl;
1651       }
1652       if (ipar == int(parResol.size() + parScale.size())) {
1653         FitParametersFile << " Cross section fit parameters:" << std::endl;
1654       }
1655       if (ipar == int(parResol.size() + parScale.size() + crossSectionHandler->parNum())) {
1656         FitParametersFile << " Background fit parameters:" << std::endl;
1657       }
1658       //       if( ipar >= int(parResol.size()+parScale.size()) && ipar < int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) && notWritted ) {
1659 
1660       //         std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]));
1661       //         std::vector<double>::const_iterator it = relativeCrossSections.begin();
1662       //         for( ; it != relativeCrossSections.end(); ++it ) {
1663       //           FitParametersFile << "  Results of the fit: parameter " << ipar << " has value "
1664       //                             << *it << "+-" << 0
1665       //                             << " + " << 0 << " - " << 0
1666       //                             << " /t/t (" << 0 << ")" << std::endl;
1667       //         }
1668 
1669       //         notWritten = false;
1670       //       }
1671       //       else {
1672       FitParametersFile << "  Results of the fit: parameter " << ipar << " has value " << pval << "+-"
1673                         << parerr[3 * ipar] << " + " << parerr[3 * ipar + 1] << " - "
1674                         << parerr[3 * ipar + 2]
1675                         // << " \t\t (" << parname[ipar] << ")"
1676                         << std::endl;
1677     }
1678     rmin.mnstat(fmin, fdem, errdef, npari, nparx, istat);  // NNBB Commented for a check!
1679     FitParametersFile << std::endl;
1680 
1681     if (minimumShapePlots_) {
1682       // Create plots of the minimum vs parameters
1683       // -----------------------------------------
1684       // Keep this after the parameters filling because it recomputes the values and it can compromise the fit results.
1685       if (somethingtodo) {
1686         std::stringstream iorderString;
1687         iorderString << iorder;
1688         std::stringstream iLoopString;
1689         iLoopString << loopCounter;
1690 
1691         for (int ipar = 0; ipar < parnumber; ipar++) {
1692           if (parfix[ipar] == 1)
1693             continue;
1694           std::cout << "plotting parameter = " << ipar + 1 << std::endl;
1695           std::stringstream iparString;
1696           iparString << ipar + 1;
1697           std::stringstream iparStringName;
1698           iparStringName << ipar;
1699           rmin.mncomd(("scan " + iparString.str()).c_str(), ierror);
1700           if (ierror == 0) {
1701             TCanvas *canvas = new TCanvas(("likelihoodCanvas_loop_" + iLoopString.str() + "_oder_" +
1702                                            iorderString.str() + "_par_" + iparStringName.str())
1703                                               .c_str(),
1704                                           ("likelihood_" + iparStringName.str()).c_str(),
1705                                           1000,
1706                                           800);
1707             canvas->cd();
1708             // arglis[0] = ipar;
1709             // rmin.mnexcm( "SCA", arglis, 0, ierror );
1710             TGraph *graph = (TGraph *)rmin.GetPlot();
1711             graph->Draw("AP");
1712             // graph->SetTitle(("parvalue["+iparStringName.str()+"]").c_str());
1713             graph->SetTitle(parname[ipar]);
1714             // graph->Write();
1715 
1716             canvas->Write();
1717           }
1718         }
1719 
1720         //       // Draw contours of the fit
1721         //       TCanvas * canvas = new TCanvas(("contourCanvas_oder_"+iorderString.str()).c_str(), "contour", 1000, 800);
1722         //       canvas->cd();
1723         //       TGraph * contourGraph = (TGraph*)rmin.Contour(4, 2, 4);
1724         //       if( (rmin.GetStatus() == 0) || (rmin.GetStatus() >= 3) ) {
1725         //         contourGraph->Draw("AP");
1726         //       }
1727         //       else {
1728         //         std::cout << "Contour graph error: status = " << rmin.GetStatus() << std::endl;
1729         //       }
1730         //       canvas->Write();
1731       }
1732     }
1733 
1734   }  // end loop on iorder
1735   FitParametersFile.close();
1736 
1737   std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters after likelihood " << std::endl;
1738   for (unsigned int ipar = 0; ipar < (unsigned int)parnumber; ipar++) {
1739     std::cout << ipar << " " << parvalue[loopCounter][ipar] << " : free = " << parfix[ipar]
1740               << "; order = " << parorder[ipar] << std::endl;
1741   }
1742 
1743   // Put back parvalue into parResol, parScale, parCrossSection, parBgr
1744   // ------------------------------------------------------------------
1745   for (int i = 0; i < (int)(parResol.size()); ++i) {
1746     parResol[i] = parvalue[loopCounter][i];
1747   }
1748   for (int i = 0; i < (int)(parScale.size()); ++i) {
1749     parScale[i] = parvalue[loopCounter][i + parResol.size()];
1750   }
1751   parCrossSection =
1752       crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size() + parScale.size()]), resfind);
1753   for (unsigned int i = 0; i < parCrossSection.size(); ++i) {
1754     //   parCrossSection[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()];
1755     std::cout << "relative cross section[" << i << "] = " << parCrossSection[i] << std::endl;
1756   }
1757   // Save only the fitted background parameters
1758   for (unsigned int i = 0; i < (parBgr.size() - parForResonanceWindows); ++i) {
1759     parBgr[i] = parvalue[loopCounter][i + parResol.size() + parScale.size() + crossSectionHandler->parNum()];
1760   }
1761 
1762   // Background rescaling from regions to resonances
1763   // -----------------------------------------------
1764   // Only if we fitted the background
1765   if (doBackgroundFit[loopCounter]) {
1766     // This rescales from regions to resonances
1767     int localMuonType = MuonType;
1768     if (MuonType > 2)
1769       localMuonType = 2;
1770     backgroundHandler->rescale(parBgr, ResMass, massWindowHalfWidth[localMuonType], MuScleFitUtils::ReducedSavedPair);
1771   }
1772 
1773   // Delete the arrays used to set some parameters
1774   delete[] Start;
1775   delete[] Step;
1776   delete[] Mini;
1777   delete[] Maxi;
1778   delete[] ind;
1779   delete[] parname;
1780 }
1781 
1782 // Likelihood function
1783 // -------------------
1784 extern "C" void likelihood(int &npar, double *grad, double &fval, double *xval, int flag) {
1785   if (MuScleFitUtils::debug > 19)
1786     std::cout << "[MuScleFitUtils-likelihood]: In likelihood function" << std::endl;
1787 
1788   const lorentzVector *recMu1;
1789   const lorentzVector *recMu2;
1790   lorentzVector corrMu1;
1791   lorentzVector corrMu2;
1792 
1793   //   if (MuScleFitUtils::debug>19) {
1794   //     int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
1795   //                           MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
1796   //     std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
1797   //     for (int ipar=0; ipar<parnumber; ipar++) {
1798   //       std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
1799   //     }
1800   //     std::cout << std::endl;
1801   //   }
1802 
1803   // Loop on the tree
1804   // ----------------
1805   double flike = 0;
1806   int evtsinlik = 0;
1807   int evtsoutlik = 0;
1808   // std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
1809   if (MuScleFitUtils::debug > 0) {
1810     std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
1811     std::cout << "ReducedSavedPair.size() = " << MuScleFitUtils::ReducedSavedPair.size() << std::endl;
1812   }
1813   // for( unsigned int nev=0; nev<MuScleFitUtils::SavedPair.size(); ++nev ) {
1814   for (unsigned int nev = 0; nev < MuScleFitUtils::ReducedSavedPair.size(); ++nev) {
1815     //     recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
1816     //     recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
1817     recMu1 = &(MuScleFitUtils::ReducedSavedPair[nev].first);
1818     recMu2 = &(MuScleFitUtils::ReducedSavedPair[nev].second);
1819 
1820     // Compute original mass
1821     // ---------------------
1822     double mass = MuScleFitUtils::invDimuonMass(*recMu1, *recMu2);
1823 
1824     // Compute weight and reference mass (from original mass)
1825     // ------------------------------------------------------
1826     double weight = MuScleFitUtils::computeWeight(mass, MuScleFitUtils::iev_);
1827     if (weight != 0.) {
1828       // Compute corrected mass (from previous biases) only if we are currently fitting the scale
1829       // ----------------------------------------------------------------------------------------
1830       if (MuScleFitUtils::doScaleFit[MuScleFitUtils::loopCounter]) {
1831         //  std::cout << "Original pt1 = " << corrMu1.Pt() << std::endl;
1832         //  std::cout << "Original pt2 = " << corrMu2.Pt() << std::endl;
1833         corrMu1 = MuScleFitUtils::applyScale(*recMu1, xval, -1);
1834         corrMu2 = MuScleFitUtils::applyScale(*recMu2, xval, 1);
1835 
1836         //         if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
1837         //           std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
1838         //           std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
1839         //         }
1840         //  std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
1841         //  std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
1842       } else {
1843         corrMu1 = *recMu1;
1844         corrMu2 = *recMu2;
1845 
1846         //         if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
1847         //           std::cout << "Not rescaled pt1 = " << corrMu1.Pt() << std::endl;
1848         //           std::cout << "Not rescaled pt2 = " << corrMu2.Pt() << std::endl;
1849         //         }
1850       }
1851       double corrMass = MuScleFitUtils::invDimuonMass(corrMu1, corrMu2);
1852       double Y = (corrMu1 + corrMu2).Rapidity();
1853       double resEta = (corrMu1 + corrMu2).Eta();
1854       if (MuScleFitUtils::debug > 19) {
1855         std::cout << "[MuScleFitUtils-likelihood]: Original/Corrected resonance mass = " << mass << " / " << corrMass
1856                   << std::endl;
1857       }
1858 
1859       // Compute mass resolution
1860       // -----------------------
1861       double massResol = MuScleFitUtils::massResolution(corrMu1, corrMu2, xval);
1862       if (MuScleFitUtils::debug > 19)
1863         std::cout << "[MuScleFitUtils-likelihood]: Resolution is " << massResol << std::endl;
1864 
1865       // Compute probability of this mass value including background modeling
1866       // --------------------------------------------------------------------
1867       if (MuScleFitUtils::debug > 1)
1868         std::cout << "calling massProb inside likelihood function" << std::endl;
1869 
1870       // double prob = MuScleFitUtils::massProb( corrMass, resEta, Y, massResol, xval );
1871       double prob = MuScleFitUtils::massProb(corrMass, resEta, Y, massResol, xval, false, corrMu1.eta(), corrMu2.eta());
1872       if (MuScleFitUtils::debug > 1)
1873         std::cout << "likelihood:massProb = " << prob << std::endl;
1874 
1875       // Compute likelihood
1876       // ------------------
1877       if (prob > 0) {
1878         // flike += log(prob*10000)*weight; // NNBB! x10000 to see if we can recover the problem of boundary
1879         flike += log(prob) * weight;
1880         evtsinlik += 1;  // NNBB test: see if likelihood per event is smarter (boundary problem)
1881       } else {
1882         if (MuScleFitUtils::debug > 0) {
1883           std::cout << "WARNING: corrMass = " << corrMass
1884                     << " outside window, this will cause a discontinuity in the likelihood. Consider increasing the "
1885                        "safety bands which are now set to 90% of the normalization window to avoid this problem"
1886                     << std::endl;
1887           std::cout << "Original mass was = " << mass << std::endl;
1888           std::cout << "WARNING: massResol = " << massResol << " outside window" << std::endl;
1889         }
1890         evtsoutlik += 1;
1891       }
1892       if (MuScleFitUtils::debug > 19)
1893         std::cout << "[MuScleFitUtils-likelihood]: Mass probability = " << prob << std::endl;
1894     }  // weight!=0
1895 
1896   }  // End of loop on tree events
1897 
1898   //   // Protection for low statistic. If the likelihood manages to throw out all the signal
1899   //   // events and stays with ~ 10 events in the resonance window it could have a better likelihood
1900   //   // because of ~ uniformly distributed events (a random combination could be good and spoil the fit).
1901   //   // We require that the number of events included in the fit does not change more than 5% in each minuit loop.
1902   //   bool lowStatPenalty = false;
1903   //   if( MuScleFitUtils::minuitLoop_ > 0 ) {
1904   //     double newEventsOutInRatio = double(evtsinlik);
1905   //     // double newEventsOutInRatio = double(evtsoutlik)/double(evtsinlik);
1906   //     double ratio = newEventsOutInRatio/MuScleFitUtils::oldEventsOutInRatio_;
1907   //     MuScleFitUtils::oldEventsOutInRatio_ = newEventsOutInRatio;
1908   //     if( ratio < 0.8 || ratio > 1.2 ) {
1909   //       std::cout << "Warning: too much change from oldEventsInLikelihood to newEventsInLikelihood, ratio is = " << ratio << std::endl;
1910   //       std::cout << "oldEventsInLikelihood = " << MuScleFitUtils::oldEventsOutInRatio_ << ", newEventsInLikelihood = " << newEventsOutInRatio << std::endl;
1911   //       lowStatPenalty = true;
1912   //     }
1913   //   }
1914 
1915   // It is a product of probabilities, we compare the sqrt_N of them. Thus N becomes a denominator of the logarithm.
1916   if (evtsinlik != 0) {
1917     if (MuScleFitUtils::normalizeLikelihoodByEventNumber_) {
1918       // && !(MuScleFitUtils::duringMinos_) ) {
1919       if (MuScleFitUtils::rminPtr_ == nullptr) {
1920         std::cout << "ERROR: rminPtr_ = " << MuScleFitUtils::rminPtr_ << ", code will crash" << std::endl;
1921       }
1922       double normalizationArg[] = {1 / double(evtsinlik)};
1923       // Reset the normalizationArg only if it changed
1924       if (MuScleFitUtils::oldNormalization_ != normalizationArg[0]) {
1925         int ierror = 0;
1926         //         if( MuScleFitUtils::likelihoodInLoop_ != 0 ) {
1927         //           // This condition is set only when minimizing. Later calls of hesse and minos will not change the value
1928         //           // This is done to avoid minos being confused by changing the UP parameter during its computation.
1929         //           MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
1930         //         }
1931         MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
1932         std::cout << "oldNormalization = " << MuScleFitUtils::oldNormalization_ << " new = " << normalizationArg[0]
1933                   << std::endl;
1934         MuScleFitUtils::oldNormalization_ = normalizationArg[0];
1935         MuScleFitUtils::normalizationChanged_ += 1;
1936       }
1937       fval = -2. * flike / double(evtsinlik);
1938       // fval = -2.*flike;
1939       //     if( lowStatPenalty ) {
1940       //       fval *= 100;
1941       //     }
1942     } else {
1943       fval = -2. * flike;
1944     }
1945   } else {
1946     std::cout << "Problem: Events in likelihood = " << evtsinlik << std::endl;
1947     fval = 999999999.;
1948   }
1949   // fval = -2.*flike;
1950   if (MuScleFitUtils::debug > 19)
1951     std::cout << "[MuScleFitUtils-likelihood]: End tree loop with likelihood value = " << fval << std::endl;
1952 
1953   //  #ifdef DEBUG
1954 
1955   //  if( MuScleFitUtils::minuitLoop_ < 10000 ) {
1956   if (MuScleFitUtils::likelihoodInLoop_ != nullptr) {
1957     ++MuScleFitUtils::minuitLoop_;
1958     MuScleFitUtils::likelihoodInLoop_->SetBinContent(MuScleFitUtils::minuitLoop_, fval);
1959   }
1960   //  }
1961   // else std::cout << "minuitLoop over 10000. Not filling histogram" << std::endl;
1962 
1963   std::cout << "MINUIT loop number " << MuScleFitUtils::minuitLoop_ << ", likelihood = " << fval << std::endl;
1964 
1965   if (MuScleFitUtils::debug > 0) {
1966     //     if( MuScleFitUtils::duringMinos_ ) {
1967     //       int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
1968     //                             MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
1969     //       std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
1970     //       for (int ipar=0; ipar<parnumber; ipar++) {
1971     //         std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
1972     //       }
1973     //       std::cout << std::endl;
1974     //       std::cout << "[MuScleFitUtils-likelihood]: likelihood value = " << fval << std::endl;
1975     //     }
1976     std::cout << "Events in likelihood = " << evtsinlik << std::endl;
1977     std::cout << "Events out likelihood = " << evtsoutlik << std::endl;
1978   }
1979 
1980   //  #endif
1981 }
1982 
1983 // Mass fitting routine
1984 // --------------------
1985 std::vector<TGraphErrors *> MuScleFitUtils::fitMass(TH2F *histo) {
1986   if (MuScleFitUtils::debug > 0)
1987     std::cout << "Fitting " << histo->GetName() << std::endl;
1988 
1989   std::vector<TGraphErrors *> results;
1990 
1991   // Results of the fit
1992   // ------------------
1993   std::vector<double> Ftop;
1994   std::vector<double> Fwidth;
1995   std::vector<double> Fmass;
1996   std::vector<double> Etop;
1997   std::vector<double> Ewidth;
1998   std::vector<double> Emass;
1999   std::vector<double> Fchi2;
2000   // X bin center and width
2001   // ----------------------
2002   std::vector<double> Xcenter;
2003   std::vector<double> Ex;
2004 
2005   // Fit with lorentzian peak
2006   // ------------------------
2007   TF1 *fitFcn = new TF1("fitFcn", lorentzianPeak, 70, 110, 3);
2008   fitFcn->SetParameters(100, 3, 91);
2009   fitFcn->SetParNames("Ftop", "Fwidth", "Fmass");
2010   fitFcn->SetLineWidth(2);
2011 
2012   // Fit slices projected along Y from bins in X
2013   // -------------------------------------------
2014   double cont_min = 20;  // Minimum number of entries
2015   Int_t binx = histo->GetXaxis()->GetNbins();
2016   // TFile *f= new TFile("prova.root", "recreate");
2017   // histo->Write();
2018   for (int i = 1; i <= binx; i++) {
2019     TH1 *histoY = histo->ProjectionY("", i, i);
2020     // histoY->Write();
2021     double cont = histoY->GetEntries();
2022     if (cont > cont_min) {
2023       histoY->Fit("fitFcn", "0", "", 70, 110);
2024       double *par = fitFcn->GetParameters();
2025       const double *err = fitFcn->GetParErrors();
2026 
2027       Ftop.push_back(par[0]);
2028       Fwidth.push_back(par[1]);
2029       Fmass.push_back(par[2]);
2030       Etop.push_back(err[0]);
2031       Ewidth.push_back(err[1]);
2032       Emass.push_back(err[2]);
2033 
2034       double chi2 = fitFcn->GetChisquare();
2035       Fchi2.push_back(chi2);
2036 
2037       double xx = histo->GetXaxis()->GetBinCenter(i);
2038       Xcenter.push_back(xx);
2039       double ex = 0;  // FIXME: you can use the bin width
2040       Ex.push_back(ex);
2041     }
2042   }
2043   // f->Close();
2044 
2045   // Put the fit results in arrays for TGraphErrors
2046   // ----------------------------------------------
2047   const int nn = Fmass.size();
2048   double *x = new double[nn];
2049   double *ym = new double[nn];
2050   double *e = new double[nn];
2051   double *eym = new double[nn];
2052   double *yw = new double[nn];
2053   double *eyw = new double[nn];
2054   double *yc = new double[nn];
2055 
2056   for (int j = 0; j < nn; j++) {
2057     x[j] = Xcenter[j];
2058     ym[j] = Fmass[j];
2059     eym[j] = Emass[j];
2060     yw[j] = Fwidth[j];
2061     eyw[j] = Ewidth[j];
2062     yc[j] = Fchi2[j];
2063     e[j] = Ex[j];
2064   }
2065 
2066   // Create TGraphErrors
2067   // -------------------
2068   TString name = histo->GetName();
2069   TGraphErrors *grM = new TGraphErrors(nn, x, ym, e, eym);
2070   grM->SetTitle(name + "_M");
2071   grM->SetName(name + "_M");
2072   TGraphErrors *grW = new TGraphErrors(nn, x, yw, e, eyw);
2073   grW->SetTitle(name + "_W");
2074   grW->SetName(name + "_W");
2075   TGraphErrors *grC = new TGraphErrors(nn, x, yc, e, e);
2076   grC->SetTitle(name + "_chi2");
2077   grC->SetName(name + "_chi2");
2078 
2079   // Cleanup
2080   // -------
2081   delete[] x;
2082   delete[] ym;
2083   delete[] eym;
2084   delete[] yw;
2085   delete[] eyw;
2086   delete[] yc;
2087   delete[] e;
2088   delete fitFcn;
2089 
2090   results.push_back(grM);
2091   results.push_back(grW);
2092   results.push_back(grC);
2093   return results;
2094 }
2095 
2096 // Resolution fitting routine
2097 // --------------------------
2098 std::vector<TGraphErrors *> MuScleFitUtils::fitReso(TH2F *histo) {
2099   std::cout << "Fitting " << histo->GetName() << std::endl;
2100   std::vector<TGraphErrors *> results;
2101 
2102   // Results from fit
2103   // ----------------
2104   std::vector<double> maxs;
2105   std::vector<double> means;
2106   std::vector<double> sigmas;
2107   std::vector<double> chi2s;
2108   std::vector<double> Emaxs;
2109   std::vector<double> Emeans;
2110   std::vector<double> Esigmas;
2111 
2112   // X bin center and width
2113   // ----------------------
2114   std::vector<double> Xcenter;
2115   std::vector<double> Ex;
2116 
2117   // Fit with a gaussian
2118   // -------------------
2119   TF1 *fitFcn = new TF1("fitFunc", Gaussian, -0.2, 0.2, 3);
2120   fitFcn->SetParameters(100, 0, 0.02);
2121   fitFcn->SetParNames("max", "mean", "sigma");
2122   fitFcn->SetLineWidth(2);
2123 
2124   // Fit slices projected along Y from bins in X
2125   // -------------------------------------------
2126   double cont_min = 20;  // Minimum number of entries
2127   Int_t binx = histo->GetXaxis()->GetNbins();
2128   for (int i = 1; i <= binx; i++) {
2129     TH1 *histoY = histo->ProjectionY("", i, i);
2130     double cont = histoY->GetEntries();
2131     if (cont > cont_min) {
2132       histoY->Fit("fitFunc", "0", "", -0.2, 0.2);
2133       double *par = fitFcn->GetParameters();
2134       const double *err = fitFcn->GetParErrors();
2135 
2136       maxs.push_back(par[0]);
2137       means.push_back(par[1]);
2138       sigmas.push_back(par[2]);
2139       Emaxs.push_back(err[0]);
2140       Emeans.push_back(err[1]);
2141       Esigmas.push_back(err[2]);
2142 
2143       double chi2 = fitFcn->GetChisquare();
2144       chi2s.push_back(chi2);
2145 
2146       double xx = histo->GetXaxis()->GetBinCenter(i);
2147       Xcenter.push_back(xx);
2148       double ex = 0;  // FIXME: you can use the bin width
2149       Ex.push_back(ex);
2150     }
2151   }
2152 
2153   // Put the fit results in arrays for TGraphErrors
2154   // ----------------------------------------------
2155   const int nn = means.size();
2156   double *x = new double[nn];
2157   double *ym = new double[nn];
2158   double *e = new double[nn];
2159   double *eym = new double[nn];
2160   double *yw = new double[nn];
2161   double *eyw = new double[nn];
2162   double *yc = new double[nn];
2163 
2164   for (int j = 0; j < nn; j++) {
2165     x[j] = Xcenter[j];
2166     ym[j] = means[j];
2167     eym[j] = Emeans[j];
2168     // yw[j]  = maxs[j];
2169     // eyw[j] = Emaxs[j];
2170     yw[j] = sigmas[j];
2171     eyw[j] = Esigmas[j];
2172     yc[j] = chi2s[j];
2173     e[j] = Ex[j];
2174   }
2175 
2176   // Create TGraphErrors
2177   // -------------------
2178   TString name = histo->GetName();
2179   TGraphErrors *grM = new TGraphErrors(nn, x, ym, e, eym);
2180   grM->SetTitle(name + "_mean");
2181   grM->SetName(name + "_mean");
2182   TGraphErrors *grW = new TGraphErrors(nn, x, yw, e, eyw);
2183   grW->SetTitle(name + "_sigma");
2184   grW->SetName(name + "_sigma");
2185   TGraphErrors *grC = new TGraphErrors(nn, x, yc, e, e);
2186   grC->SetTitle(name + "_chi2");
2187   grC->SetName(name + "_chi2");
2188 
2189   // Cleanup
2190   // -------
2191   delete[] x;
2192   delete[] ym;
2193   delete[] eym;
2194   delete[] yw;
2195   delete[] eyw;
2196   delete[] yc;
2197   delete[] e;
2198   delete fitFcn;
2199 
2200   results.push_back(grM);
2201   results.push_back(grW);
2202   results.push_back(grC);
2203   return results;
2204 }
2205 
2206 // Mass probability for likelihood computation - no-background version (not used anymore)
2207 // --------------------------------------------------------------------------------------
2208 double MuScleFitUtils::massProb(const double &mass, const double &rapidity, const int ires, const double &massResol) {
2209   // This routine computes the likelihood that a given measured mass "measMass" is
2210   // the result of resonance #ires if the resolution expected for the two muons is massResol
2211   // ---------------------------------------------------------------------------------------
2212 
2213   double P = 0.;
2214 
2215   // Return Lorentz value for zero resolution cases (like MC)
2216   // --------------------------------------------------------
2217   if (massResol == 0.) {
2218     if (debug > 9)
2219       std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2220                 << massResol << " : used Lorentz P-value" << std::endl;
2221     return (0.5 * ResGamma[ires] / TMath::Pi()) /
2222            ((mass - ResMass[ires]) * (mass - ResMass[ires]) + .25 * ResGamma[ires] * ResGamma[ires]);
2223   }
2224 
2225   // NB defined as below, P is not a "probability" but a likelihood that we observe
2226   // a dimuon mass "mass", given mRef, gamma, and massResol. It is what we need for the
2227   // fit which finds the best resolution parameters, though. A definition which is
2228   // more properly a probability is given below (in massProb2()), where the result
2229   // cannot be used to fit resolution parameters because the fit would always prefer
2230   // to set the res parameters to the minimum possible value (best resolution),
2231   // to have a probability close to one of observing any mass.
2232   // -------------------------------------------------------------------------------
2233   // NNBB the following two lines have been replaced with the six following them,
2234   // which provide an improvement of a factor 9 in speed of execution at a
2235   // negligible price in precision.
2236   // ----------------------------------------------------------------------------
2237   // GL->SetParameters(gamma,mRef,mass,massResol);
2238   // P = GL->Integral(mass-5*massResol, mass+5*massResol);
2239 
2240   Int_t np = 100;
2241   double *x = new double[np];
2242   double *w = new double[np];
2243   GL->SetParameters(ResGamma[ires], ResMass[ires], mass, massResol);
2244   GL->CalcGaussLegendreSamplingPoints(np, x, w, 0.1e-15);
2245   P = GL->IntegralFast(np, x, w, ResMass[ires] - 10 * ResGamma[ires], ResMass[ires] + 10 * ResGamma[ires]);
2246   delete[] x;
2247   delete[] w;
2248 
2249   // If we are too far away we set P to epsilon and forget about this event
2250   // ----------------------------------------------------------------------
2251   if (P < 1.0e-12) {
2252     P = 1.0e-12;
2253     if (debug > 9)
2254       std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2255                 << massResol << ": used epsilon" << std::endl;
2256     return P;
2257   }
2258 
2259   if (debug > 9)
2260     std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2261               << massResol << " " << P << std::endl;
2262   return P;
2263 }
2264 
2265 std::pair<lorentzVector, lorentzVector> MuScleFitUtils::findSimMuFromRes(
2266     const edm::Handle<edm::HepMCProduct> &evtMC, const edm::Handle<edm::SimTrackContainer> &simTracks) {
2267   //Loop on simulated tracks
2268   std::pair<lorentzVector, lorentzVector> simMuFromRes;
2269   for (edm::SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
2270     //Chose muons
2271     if (std::abs((*simTrack).type()) == 13) {
2272       //If tracks from IP than find mother
2273       if ((*simTrack).genpartIndex() > 0) {
2274         HepMC::GenParticle *gp = evtMC->GetEvent()->barcode_to_particle((*simTrack).genpartIndex());
2275         if (gp != nullptr) {
2276           for (HepMC::GenVertex::particle_iterator mother = gp->production_vertex()->particles_begin(HepMC::ancestors);
2277                mother != gp->production_vertex()->particles_end(HepMC::ancestors);
2278                ++mother) {
2279             bool fromRes = false;
2280             unsigned int motherPdgId = (*mother)->pdg_id();
2281             for (int ires = 0; ires < 6; ++ires) {
2282               if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2283                 fromRes = true;
2284             }
2285             if (fromRes) {
2286               if (gp->pdg_id() == 13)
2287                 simMuFromRes.first = lorentzVector(simTrack->momentum().px(),
2288                                                    simTrack->momentum().py(),
2289                                                    simTrack->momentum().pz(),
2290                                                    simTrack->momentum().e());
2291               else
2292                 simMuFromRes.second = lorentzVector(simTrack->momentum().px(),
2293                                                     simTrack->momentum().py(),
2294                                                     simTrack->momentum().pz(),
2295                                                     simTrack->momentum().e());
2296             }
2297           }
2298         }
2299         // else LogDebug("MuScleFitUtils") << "WARNING: no matching genParticle found for simTrack" << std::endl;
2300       }
2301     }
2302   }
2303   return simMuFromRes;
2304 }
2305 
2306 std::pair<lorentzVector, lorentzVector> MuScleFitUtils::findGenMuFromRes(const edm::HepMCProduct *evtMC) {
2307   const HepMC::GenEvent *Evt = evtMC->GetEvent();
2308   std::pair<lorentzVector, lorentzVector> muFromRes;
2309   //Loop on generated particles
2310   for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); part++) {
2311     if (std::abs((*part)->pdg_id()) == 13 && (*part)->status() == 1) {
2312       bool fromRes = false;
2313       for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);
2314            mother != (*part)->production_vertex()->particles_end(HepMC::ancestors);
2315            ++mother) {
2316         unsigned int motherPdgId = (*mother)->pdg_id();
2317 
2318         // For sherpa the resonance is not saved. The muons from the resonance can be identified
2319         // by having as mother a muon of status 3.
2320         if (sherpa_) {
2321           if (motherPdgId == 13 && (*mother)->status() == 3)
2322             fromRes = true;
2323         } else {
2324           for (int ires = 0; ires < 6; ++ires) {
2325             if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2326               fromRes = true;
2327           }
2328         }
2329       }
2330       if (fromRes) {
2331         if ((*part)->pdg_id() == 13)
2332           //   muFromRes.first = (*part)->momentum();
2333           muFromRes.first = (lorentzVector(
2334               (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
2335         else
2336           muFromRes.second = (lorentzVector(
2337               (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
2338       }
2339     }
2340   }
2341   return muFromRes;
2342 }
2343 
2344 std::pair<lorentzVector, lorentzVector> MuScleFitUtils::findGenMuFromRes(
2345     const reco::GenParticleCollection *genParticles) {
2346   std::pair<lorentzVector, lorentzVector> muFromRes;
2347 
2348   //Loop on generated particles
2349   if (debug > 0)
2350     std::cout << "Starting loop on " << genParticles->size() << " genParticles" << std::endl;
2351   for (reco::GenParticleCollection::const_iterator part = genParticles->begin(); part != genParticles->end(); ++part) {
2352     if (std::abs(part->pdgId()) == 13 && part->status() == 1) {
2353       bool fromRes = false;
2354       unsigned int motherPdgId = part->mother()->pdgId();
2355       if (debug > 0) {
2356         std::cout << "Found a muon with mother: " << motherPdgId << std::endl;
2357       }
2358       for (int ires = 0; ires < 6; ++ires) {
2359         if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2360           fromRes = true;
2361       }
2362       if (fromRes) {
2363         if (part->pdgId() == 13) {
2364           muFromRes.first = part->p4();
2365           if (debug > 0)
2366             std::cout << "Found a genMuon + : " << muFromRes.first << std::endl;
2367           //      muFromRes.first = (lorentzVector(part->p4().px(),part->p4().py(),
2368           //                       part->p4().pz(),part->p4().e()));
2369         } else {
2370           muFromRes.second = part->p4();
2371           if (debug > 0)
2372             std::cout << "Found a genMuon - : " << muFromRes.second << std::endl;
2373           //      muFromRes.second = (lorentzVector(part->p4().px(),part->p4().py(),
2374           //                        part->p4().pz(),part->p4().e()));
2375         }
2376       }
2377     }
2378   }
2379   return muFromRes;
2380 }