CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Classes | Typedefs | Functions
MuScleFitUtils.h File Reference
#include <CLHEP/Vector/LorentzVector.h>
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "TGraphErrors.h"
#include "TH2F.h"
#include "TMinuit.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/CrossSectionHandler.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundHandler.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/GenMuonPair.h"
#include <vector>

Go to the source code of this file.

Classes

class  biasFunctionBase< T >
 
struct  MuScleFitUtils::byPt
 
struct  MuScleFitUtils::massResolComponentsStruct
 
class  MuScleFitUtils
 
class  resolutionFunctionBase< T >
 
class  scaleFunctionBase< T >
 

Typedefs

typedef
reco::Particle::LorentzVector 
lorentzVector
 

Functions

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

Typedef Documentation

Definition at line 46 of file MuScleFitUtils.h.

Function Documentation

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

Definition at line 1703 of file MuScleFitUtils.cc.

References MuScleFitUtils::applyScale(), MuScleFitUtils::computeWeight(), gather_cfg::cout, MuScleFitUtils::debug, MuScleFitUtils::doScaleFit, reco::tau::disc::Eta(), MuScleFitUtils::iev_, MuScleFitUtils::invDimuonMass(), MuScleFitUtils::likelihoodInLoop_, create_public_lumi_plots::log, MuScleFitUtils::loopCounter, scaleCards::mass, MuScleFitUtils::massProb(), MuScleFitUtils::massResolution(), MuScleFitUtils::minuitLoop_, MuScleFitUtils::normalizationChanged_, MuScleFitUtils::normalizeLikelihoodByEventNumber_, MuScleFitUtils::oldNormalization_, MuScleFitUtils::ReducedSavedPair, MuScleFitUtils::rminPtr_, MuScleFitUtils::SavedPair, and CommonMethods::weight().

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

1703  {
1704 
1705  if (MuScleFitUtils::debug>19) std::cout << "[MuScleFitUtils-likelihood]: In likelihood function" << std::endl;
1706 
1707  const lorentzVector * recMu1;
1708  const lorentzVector * recMu2;
1709  lorentzVector corrMu1;
1710  lorentzVector corrMu2;
1711 
1712  // if (MuScleFitUtils::debug>19) {
1713  // int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
1714  // MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
1715  // std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
1716  // for (int ipar=0; ipar<parnumber; ipar++) {
1717  // std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
1718  // }
1719  // std::cout << std::endl;
1720  // }
1721 
1722  // Loop on the tree
1723  // ----------------
1724  double flike = 0;
1725  int evtsinlik = 0;
1726  int evtsoutlik = 0;
1727  // std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
1728  if( MuScleFitUtils::debug>0 ) {
1729  std::cout << "SavedPair.size() = " << MuScleFitUtils::SavedPair.size() << std::endl;
1730  std::cout << "ReducedSavedPair.size() = " << MuScleFitUtils::ReducedSavedPair.size() << std::endl;
1731  }
1732  // for( unsigned int nev=0; nev<MuScleFitUtils::SavedPair.size(); ++nev ) {
1733  for( unsigned int nev=0; nev<MuScleFitUtils::ReducedSavedPair.size(); ++nev ) {
1734 
1735  // recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
1736  // recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
1737  recMu1 = &(MuScleFitUtils::ReducedSavedPair[nev].first);
1738  recMu2 = &(MuScleFitUtils::ReducedSavedPair[nev].second);
1739 
1740  // Compute original mass
1741  // ---------------------
1742  double mass = MuScleFitUtils::invDimuonMass( *recMu1, *recMu2 );
1743 
1744  // Compute weight and reference mass (from original mass)
1745  // ------------------------------------------------------
1747  if( weight!=0. ) {
1748  // Compute corrected mass (from previous biases) only if we are currently fitting the scale
1749  // ----------------------------------------------------------------------------------------
1751 // std::cout << "Original pt1 = " << corrMu1.Pt() << std::endl;
1752 // std::cout << "Original pt2 = " << corrMu2.Pt() << std::endl;
1753  corrMu1 = MuScleFitUtils::applyScale(*recMu1, xval, -1);
1754  corrMu2 = MuScleFitUtils::applyScale(*recMu2, xval, 1);
1755 
1756 // if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
1757 // std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
1758 // std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
1759 // }
1760 // std::cout << "Rescaled pt1 = " << corrMu1.Pt() << std::endl;
1761 // std::cout << "Rescaled pt2 = " << corrMu2.Pt() << std::endl;
1762  }
1763  else {
1764  corrMu1 = *recMu1;
1765  corrMu2 = *recMu2;
1766 
1767 // if( (corrMu1.Pt() != corrMu1.Pt()) || (corrMu2.Pt() != corrMu2.Pt()) ) {
1768 // std::cout << "Not rescaled pt1 = " << corrMu1.Pt() << std::endl;
1769 // std::cout << "Not rescaled pt2 = " << corrMu2.Pt() << std::endl;
1770 // }
1771  }
1772  double corrMass = MuScleFitUtils::invDimuonMass(corrMu1, corrMu2);
1773  double Y = (corrMu1+corrMu2).Rapidity();
1774  double resEta = (corrMu1+corrMu2).Eta();
1775  if( MuScleFitUtils::debug>19 ) {
1776  std::cout << "[MuScleFitUtils-likelihood]: Original/Corrected resonance mass = " << mass
1777  << " / " << corrMass << std::endl;
1778  }
1779 
1780  // Compute mass resolution
1781  // -----------------------
1782  double massResol = MuScleFitUtils::massResolution(corrMu1, corrMu2, xval);
1783  if (MuScleFitUtils::debug>19)
1784  std::cout << "[MuScleFitUtils-likelihood]: Resolution is " << massResol << std::endl;
1785 
1786  // Compute probability of this mass value including background modeling
1787  // --------------------------------------------------------------------
1788  if (MuScleFitUtils::debug>1) std::cout << "calling massProb inside likelihood function" << std::endl;
1789 
1790  // double prob = MuScleFitUtils::massProb( corrMass, resEta, Y, massResol, xval );
1791  double prob = MuScleFitUtils::massProb( corrMass, resEta, Y, massResol, xval, false, corrMu1.eta(), corrMu2.eta() );
1792  if (MuScleFitUtils::debug>1) std::cout << "likelihood:massProb = " << prob << std::endl;
1793 
1794  // Compute likelihood
1795  // ------------------
1796  if( prob>0 ) {
1797  // flike += log(prob*10000)*weight; // NNBB! x10000 to see if we can recover the problem of boundary
1798  flike += log(prob)*weight;
1799  evtsinlik += 1; // NNBB test: see if likelihood per event is smarter (boundary problem)
1800  } else {
1801  if( MuScleFitUtils::debug > 0 ) {
1802  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;
1803  std::cout << "Original mass was = " << mass << std::endl;
1804  std::cout << "WARNING: massResol = " << massResol << " outside window" << std::endl;
1805  }
1806  evtsoutlik += 1;
1807  }
1808  if (MuScleFitUtils::debug>19)
1809  std::cout << "[MuScleFitUtils-likelihood]: Mass probability = " << prob << std::endl;
1810  } // weight!=0
1811 
1812  } // End of loop on tree events
1813 
1814 // // Protection for low statistic. If the likelihood manages to throw out all the signal
1815 // // events and stays with ~ 10 events in the resonance window it could have a better likelihood
1816 // // because of ~ uniformly distributed events (a random combination could be good and spoil the fit).
1817 // // We require that the number of events included in the fit does not change more than 5% in each minuit loop.
1818 // bool lowStatPenalty = false;
1819 // if( MuScleFitUtils::minuitLoop_ > 0 ) {
1820 // double newEventsOutInRatio = double(evtsinlik);
1821 // // double newEventsOutInRatio = double(evtsoutlik)/double(evtsinlik);
1822 // double ratio = newEventsOutInRatio/MuScleFitUtils::oldEventsOutInRatio_;
1823 // MuScleFitUtils::oldEventsOutInRatio_ = newEventsOutInRatio;
1824 // if( ratio < 0.8 || ratio > 1.2 ) {
1825 // std::cout << "Warning: too much change from oldEventsInLikelihood to newEventsInLikelihood, ratio is = " << ratio << std::endl;
1826 // std::cout << "oldEventsInLikelihood = " << MuScleFitUtils::oldEventsOutInRatio_ << ", newEventsInLikelihood = " << newEventsOutInRatio << std::endl;
1827 // lowStatPenalty = true;
1828 // }
1829 // }
1830 
1831  // It is a product of probabilities, we compare the sqrt_N of them. Thus N becomes a denominator of the logarithm.
1832  if( evtsinlik != 0 ) {
1833 
1835  // && !(MuScleFitUtils::duringMinos_) ) {
1836  if( MuScleFitUtils::rminPtr_ == 0 ) {
1837  std::cout << "ERROR: rminPtr_ = " << MuScleFitUtils::rminPtr_ << ", code will crash" << std::endl;
1838  }
1839  double normalizationArg[] = {1/double(evtsinlik)};
1840  // Reset the normalizationArg only if it changed
1841  if( MuScleFitUtils::oldNormalization_ != normalizationArg[0] ) {
1842  int ierror = 0;
1843 // if( MuScleFitUtils::likelihoodInLoop_ != 0 ) {
1844 // // This condition is set only when minimizing. Later calls of hesse and minos will not change the value
1845 // // This is done to avoid minos being confused by changing the UP parameter during its computation.
1846 // MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
1847 // }
1848  MuScleFitUtils::rminPtr_->mnexcm("SET ERR", normalizationArg, 1, ierror);
1849  std::cout << "oldNormalization = " << MuScleFitUtils::oldNormalization_ << " new = " << normalizationArg[0] << std::endl;
1850  MuScleFitUtils::oldNormalization_ = normalizationArg[0];
1852  }
1853  fval = -2.*flike/double(evtsinlik);
1854  // fval = -2.*flike;
1855  // if( lowStatPenalty ) {
1856  // fval *= 100;
1857  // }
1858  }
1859  else {
1860  fval = -2.*flike;
1861  }
1862  }
1863  else {
1864  std::cout << "Problem: Events in likelihood = " << evtsinlik << std::endl;
1865  fval = 999999999.;
1866  }
1867  // fval = -2.*flike;
1868  if (MuScleFitUtils::debug>19)
1869  std::cout << "[MuScleFitUtils-likelihood]: End tree loop with likelihood value = " << fval << std::endl;
1870 
1871 // #ifdef DEBUG
1872 
1873 // if( MuScleFitUtils::minuitLoop_ < 10000 ) {
1877  }
1878  // }
1879  // else std::cout << "minuitLoop over 10000. Not filling histogram" << std::endl;
1880 
1881  std::cout<<"MINUIT loop number "<<MuScleFitUtils::minuitLoop_<<", likelihood = "<<fval<<std::endl;
1882 
1883  if( MuScleFitUtils::debug > 0 ) {
1884  // if( MuScleFitUtils::duringMinos_ ) {
1885  // int parnumber = (int)(MuScleFitUtils::parResol.size()+MuScleFitUtils::parScale.size()+
1886  // MuScleFitUtils::parCrossSection.size()+MuScleFitUtils::parBgr.size());
1887  // std::cout << "[MuScleFitUtils-likelihood]: Looping on tree with ";
1888  // for (int ipar=0; ipar<parnumber; ipar++) {
1889  // std::cout << "Parameter #" << ipar << " with value " << xval[ipar] << " ";
1890  // }
1891  // std::cout << std::endl;
1892  // std::cout << "[MuScleFitUtils-likelihood]: likelihood value = " << fval << std::endl;
1893  // }
1894  std::cout << "Events in likelihood = " << evtsinlik << std::endl;
1895  std::cout << "Events out likelihood = " << evtsoutlik << std::endl;
1896  }
1897 
1898 // #endif
1899 }
static std::vector< int > doScaleFit
static unsigned int loopCounter
static int debug
static unsigned int normalizationChanged_
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static std::vector< std::pair< lorentzVector, lorentzVector > > ReducedSavedPair
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
static int minuitLoop_
static double computeWeight(const double &mass, const int iev, const bool doUseBkgrWindow=false)
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
static double invDimuonMass(const lorentzVector &mu1, const lorentzVector &mu2)
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)
static TMinuit * rminPtr_
static TH1D * likelihoodInLoop_
static double oldNormalization_
static int iev_
static bool normalizeLikelihoodByEventNumber_
tuple mass
Definition: scaleCards.py:27
tuple cout
Definition: gather_cfg.py:121