#include "MuScleFitUtils.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "TString.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TF1.h"
#include "TF2.h"
#include <iostream>
#include <fstream>
#include <memory>
#include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"

Functions
Double_t	Gaussian (Double_t x, Double_t par)

void	likelihood (int &npar, double grad, double &fval, double xval, int flag)

Double_t	lorentzianPeak (Double_t x, Double_t par)

Variables
double	f [11][100]

double	g [11][100]

TF1 *	GL

TF2 *	GL2

double	isum

double	mzsum

Function Documentation

◆ Gaussian()

Double_t Gaussian	(	Double_t *	x,
		Double_t *	par
	)

Definition at line 63 of file MuScleFitUtils.cc.

References JetChargeProducer_cfi::exp, and x.

Referenced by MuScleFitUtils::fitReso().

                                               {
   return par[0] * exp(-0.5 * ((x[0] - par[1]) / par[2]) * ((x[0] - par[1]) / par[2]));
 }

◆ likelihood()

void likelihood	(	int &	npar,
		double *	grad,
		double &	fval,
		double *	xval,
		int	flag
	)

Definition at line 1784 of file MuScleFitUtils.cc.

References MuScleFitUtils::applyScale(), MuScleFitUtils::computeWeight(), gather_cfg::cout, MuScleFitUtils::debug, MuScleFitUtils::doScaleFit, MuScleFitUtils::iev_, MuScleFitUtils::invDimuonMass(), MuScleFitUtils::likelihoodInLoop_, CrabHelper::log, MuScleFitUtils::loopCounter, EgHLTOffHistBins_cfi::mass, MuScleFitUtils::massProb(), MuScleFitUtils::massResolution(), MuScleFitUtils::minuitLoop_, makeMEIFBenchmarkPlots::nev, MuScleFitUtils::normalizationChanged_, MuScleFitUtils::normalizeLikelihoodByEventNumber_, MuScleFitUtils::oldNormalization_, TtFullHadEvtBuilder_cfi::prob, MuScleFitUtils::ReducedSavedPair, MuScleFitUtils::rminPtr_, MuScleFitUtils::SavedPair, mps_merge::weight, makeMEIFBenchmarkPlots::xval, and beamSpotPI::Y.

Referenced by MuScleFitUtils::minimizeLikelihood(), and PFNuclearProducer::produce().

                                                                                           {
   if (MuScleFitUtils::debug > 19)
     std::cout << "[MuScleFitUtils-likelihood]: In likelihood function" << std::endl;
 
   const lorentzVector *recMu1;
   const lorentzVector *recMu2;
   lorentzVector corrMu1;
   lorentzVector corrMu2;
 
   //   if (MuScleFitUtils::debug>19) {
   //     int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
   //                           MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
   //     std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
   //     for (int ipar=0; ipar<parnumber; ipar++) {
   //       std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
   //     }
   //     std::cout << std::endl;
   //   }
 
   // Loop on the tree
   // ----------------
   double flike = 0;
   int evtsinlik = 0;
   int evtsoutlik = 0;
   // std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
   if (MuScleFitUtils::debug > 0) {
     std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
     std::cout << "ReducedSavedPair.size() = " << MuScleFitUtils::ReducedSavedPair.size() << std::endl;
   }
   // for( unsigned int nev=0; nev<MuScleFitUtils::SavedPair.size(); ++nev ) {
   for (unsigned int nev = 0; nev < MuScleFitUtils::ReducedSavedPair.size(); ++nev) {
     //     recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
     //     recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
     recMu1 = &(MuScleFitUtils::ReducedSavedPair[nev].first);
     recMu2 = &(MuScleFitUtils::ReducedSavedPair[nev].second);
 
     // Compute original mass
     // ---------------------
     double mass = MuScleFitUtils::invDimuonMass(*recMu1, *recMu2);
 
     // Compute weight and reference mass (from original mass)
     // ------------------------------------------------------
     double weight = MuScleFitUtils::computeWeight(mass, MuScleFitUtils::iev_);
     if (weight != 0.) {
       // Compute corrected mass (from previous biases) only if we are currently fitting the scale
       // ----------------------------------------------------------------------------------------
       if (MuScleFitUtils::doScaleFit[MuScleFitUtils::loopCounter]) {
         //  std::cout << "Original pt1 = " << corrMu1.Pt() << std::endl;
         //  std::cout << "Original pt2 = " << corrMu2.Pt() << std::endl;
         corrMu1 = MuScleFitUtils::applyScale(*recMu1, xval, -1);
         corrMu2 = MuScleFitUtils::applyScale(*recMu2, xval, 1);
 
         //         if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
         //           std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
         //           std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
         //         }
         //  std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
         //  std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
       } else {
         corrMu1 = *recMu1;
         corrMu2 = *recMu2;
 
         //         if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
         //           std::cout << "Not rescaled pt1 = " << corrMu1.Pt() << std::endl;
         //           std::cout << "Not rescaled pt2 = " << corrMu2.Pt() << std::endl;
         //         }
       }
       double corrMass = MuScleFitUtils::invDimuonMass(corrMu1, corrMu2);
       double Y = (corrMu1 + corrMu2).Rapidity();
       double resEta = (corrMu1 + corrMu2).Eta();
       if (MuScleFitUtils::debug > 19) {
         std::cout << "[MuScleFitUtils-likelihood]: Original/Corrected resonance mass = " << mass << " / " << corrMass
                   << std::endl;
       }
 
       // Compute mass resolution
       // -----------------------
       double massResol = MuScleFitUtils::massResolution(corrMu1, corrMu2, xval);
       if (MuScleFitUtils::debug > 19)
         std::cout << "[MuScleFitUtils-likelihood]: Resolution is " << massResol << std::endl;
 
       // Compute probability of this mass value including background modeling
       // --------------------------------------------------------------------
       if (MuScleFitUtils::debug > 1)
         std::cout << "calling massProb inside likelihood function" << std::endl;
 
       // double prob = MuScleFitUtils::massProb( corrMass, resEta, Y, massResol, xval );
       double prob = MuScleFitUtils::massProb(corrMass, resEta, Y, massResol, xval, false, corrMu1.eta(), corrMu2.eta());
       if (MuScleFitUtils::debug > 1)
         std::cout << "likelihood:massProb = " << prob << std::endl;
 
       // Compute likelihood
       // ------------------
       if (prob > 0) {
         // flike += log(prob*10000)*weight; // NNBB! x10000 to see if we can recover the problem of boundary
         flike += log(prob) * weight;
         evtsinlik += 1;  // NNBB test: see if likelihood per event is smarter (boundary problem)
       } else {
         if (MuScleFitUtils::debug > 0) {
           std::cout << "WARNING: corrMass = " << corrMass
                     << " outside window, this will cause a discontinuity in the likelihood. Consider increasing the "
                        "safety bands which are now set to 90% of the normalization window to avoid this problem"
                     << std::endl;
           std::cout << "Original mass was = " << mass << std::endl;
           std::cout << "WARNING: massResol = " << massResol << " outside window" << std::endl;
         }
         evtsoutlik += 1;
       }
       if (MuScleFitUtils::debug > 19)
         std::cout << "[MuScleFitUtils-likelihood]: Mass probability = " << prob << std::endl;
     }  // weight!=0
 
   }  // End of loop on tree events
 
   //   // Protection for low statistic. If the likelihood manages to throw out all the signal
   //   // events and stays with ~ 10 events in the resonance window it could have a better likelihood
   //   // because of ~ uniformly distributed events (a random combination could be good and spoil the fit).
   //   // We require that the number of events included in the fit does not change more than 5% in each minuit loop.
   //   bool lowStatPenalty = false;
   //   if( MuScleFitUtils::minuitLoop_ > 0 ) {
   //     double newEventsOutInRatio = double(evtsinlik);
   //     // double newEventsOutInRatio = double(evtsoutlik)/double(evtsinlik);
   //     double ratio = newEventsOutInRatio/MuScleFitUtils::oldEventsOutInRatio_;
   //     MuScleFitUtils::oldEventsOutInRatio_ = newEventsOutInRatio;
   //     if( ratio < 0.8 || ratio > 1.2 ) {
   //       std::cout << "Warning: too much change from oldEventsInLikelihood to newEventsInLikelihood, ratio is = " << ratio << std::endl;
   //       std::cout << "oldEventsInLikelihood = " << MuScleFitUtils::oldEventsOutInRatio_ << ", newEventsInLikelihood = " << newEventsOutInRatio << std::endl;
   //       lowStatPenalty = true;
   //     }
   //   }
 
   // It is a product of probabilities, we compare the sqrt_N of them. Thus N becomes a denominator of the logarithm.
   if (evtsinlik != 0) {
     if (MuScleFitUtils::normalizeLikelihoodByEventNumber_) {
       // && !(MuScleFitUtils::duringMinos_) ) {
       if (MuScleFitUtils::rminPtr_ == nullptr) {
         std::cout << "ERROR: rminPtr_ = " << MuScleFitUtils::rminPtr_ << ", code will crash" << std::endl;
       }
       double normalizationArg[] = {1 / double(evtsinlik)};
       // Reset the normalizationArg only if it changed
       if (MuScleFitUtils::oldNormalization_ != normalizationArg[0]) {
         int ierror = 0;
         //         if( MuScleFitUtils::likelihoodInLoop_ != 0 ) {
         //           // This condition is set only when minimizing. Later calls of hesse and minos will not change the value
         //           // This is done to avoid minos being confused by changing the UP parameter during its computation.
         //           MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
         //         }
         MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
         std::cout << "oldNormalization = " << MuScleFitUtils::oldNormalization_ << " new = " << normalizationArg[0]
                   << std::endl;
         MuScleFitUtils::oldNormalization_ = normalizationArg[0];
         MuScleFitUtils::normalizationChanged_ += 1;
       }
       fval = -2. * flike / double(evtsinlik);
       // fval = -2.*flike;
       //     if( lowStatPenalty ) {
       //       fval *= 100;
       //     }
     } else {
       fval = -2. * flike;
     }
   } else {
     std::cout << "Problem: Events in likelihood = " << evtsinlik << std::endl;
     fval = 999999999.;
   }
   // fval = -2.*flike;
   if (MuScleFitUtils::debug > 19)
     std::cout << "[MuScleFitUtils-likelihood]: End tree loop with likelihood value = " << fval << std::endl;
 
   //  #ifdef DEBUG
 
   //  if( MuScleFitUtils::minuitLoop_ < 10000 ) {
   if (MuScleFitUtils::likelihoodInLoop_ != nullptr) {
     ++MuScleFitUtils::minuitLoop_;
     MuScleFitUtils::likelihoodInLoop_->SetBinContent(MuScleFitUtils::minuitLoop_, fval);
   }
   //  }
   // else std::cout << "minuitLoop over 10000. Not filling histogram" << std::endl;
 
   std::cout << "MINUIT loop number " << MuScleFitUtils::minuitLoop_ << ", likelihood = " << fval << std::endl;
 
   if (MuScleFitUtils::debug > 0) {
     //     if( MuScleFitUtils::duringMinos_ ) {
     //       int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
     //                             MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
     //       std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
     //       for (int ipar=0; ipar<parnumber; ipar++) {
     //         std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
     //       }
     //       std::cout << std::endl;
     //       std::cout << "[MuScleFitUtils-likelihood]: likelihood value = " << fval << std::endl;
     //     }
     std::cout << "Events in likelihood = " << evtsinlik << std::endl;
     std::cout << "Events out likelihood = " << evtsoutlik << std::endl;
   }
 
   //  #endif
 }

◆ lorentzianPeak()

Double_t lorentzianPeak	(	Double_t *	x,
		Double_t *	par
	)

See header file for a class description

Author: S. Bolognesi - INFN Torino / T. Dorigo, M. De Mattia - INFN Padova Revised S. Casasso, E. Migliore - UniTo & INFN Torino

Definition at line 56 of file MuScleFitUtils.cc.

References MillePedeFileConverter_cfg::e, METSkim_cff::Max, Pi, and x.

Referenced by MuScleFitUtils::fitMass().

                                                     {
   return (0.5 * par[0] * par[1] / TMath::Pi()) /
          TMath::Max(1.e-10, (x[0] - par[2]) * (x[0] - par[2]) + .25 * par[1] * par[1]);
 }

Variable Documentation

◆ f

double f[11][100]

Definition at line 78 of file MuScleFitUtils.cc.

◆ g

double g[11][100]

Definition at line 79 of file MuScleFitUtils.cc.

◆ GL

TF1* GL

Initial value:

= new TF1(

"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)

Definition at line 83 of file MuScleFitUtils.cc.

Referenced by MuScleFitUtils::massProb(), and MuScleFitBase::readProbabilityDistributionsFromFile().

◆ GL2

TF2* GL2

Initial value:

= new TF2("GL2",
                   "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))",
                   0,
                   200,
                   0,
                   200)

Definition at line 86 of file MuScleFitUtils.cc.

◆ isum

double isum

Definition at line 77 of file MuScleFitUtils.cc.

◆ mzsum

double mzsum

Definition at line 76 of file MuScleFitUtils.cc.

Functions

Variables