CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Functions | Variables
MuScleFitUtils.cc File Reference
#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"

Go to the source code of this file.

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

Double_t Gaussian ( Double_t *  x,
Double_t *  par 
)

Definition at line 63 of file MuScleFitUtils.cc.

References funct::exp().

Referenced by MuScleFitUtils::fitReso().

63  {
64  return par[0] * exp(-0.5 * ((x[0] - par[1]) / par[2]) * ((x[0] - par[1]) / par[2]));
65 }
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
gpuClustering::x
uint16_t const *__restrict__ x
Definition: gpuClustering.h:43
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_, log, MuScleFitUtils::loopCounter, ResonanceBuilder::mass, MuScleFitUtils::massProb(), MuScleFitUtils::massResolution(), MuScleFitUtils::minuitLoop_, makeMEIFBenchmarkPlots::nev, MuScleFitUtils::normalizationChanged_, MuScleFitUtils::normalizeLikelihoodByEventNumber_, MuScleFitUtils::oldNormalization_, MuScleFitUtils::ReducedSavedPair, MuScleFitUtils::rminPtr_, MuScleFitUtils::SavedPair, histoStyle::weight, and BeamSpotPI::Y.

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

1784  {
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 }
MuScleFitUtils::doScaleFit
static std::vector< int > doScaleFit
Definition: MuScleFitUtils.h:179
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
MuScleFitUtils::loopCounter
static unsigned int loopCounter
Definition: MuScleFitUtils.h:145
MuScleFitUtils::debug
static int debug
Definition: MuScleFitUtils.h:130
MuScleFitUtils::normalizationChanged_
static unsigned int normalizationChanged_
Definition: MuScleFitUtils.h:245
MuScleFitUtils::massProb
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
Definition: MuScleFitUtils.cc:2208
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
MuScleFitUtils::ReducedSavedPair
static std::vector< std::pair< lorentzVector, lorentzVector > > ReducedSavedPair
Definition: MuScleFitUtils.h:232
MuScleFitUtils::massResolution
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
MuScleFitUtils::minuitLoop_
static int minuitLoop_
Definition: MuScleFitUtils.h:183
MuScleFitUtils::computeWeight
static double computeWeight(const double &mass, const int iev, const bool doUseBkgrWindow=false)
Definition: MuScleFitUtils.cc:1163
MuScleFitUtils::SavedPair
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
Definition: MuScleFitUtils.h:231
MuScleFitUtils::invDimuonMass
static double invDimuonMass(const lorentzVector &mu1, const lorentzVector &mu2)
Definition: MuScleFitUtils.cc:518
MuScleFitUtils::applyScale
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)
Definition: MuScleFitUtils.cc:465
MuScleFitUtils::rminPtr_
static TMinuit * rminPtr_
Definition: MuScleFitUtils.h:242
MuScleFitUtils::likelihoodInLoop_
static TH1D * likelihoodInLoop_
Definition: MuScleFitUtils.h:184
MuScleFitUtils::oldNormalization_
static double oldNormalization_
Definition: MuScleFitUtils.h:244
MuScleFitUtils::iev_
static int iev_
Definition: MuScleFitUtils.h:253
MuScleFitUtils::normalizeLikelihoodByEventNumber_
static bool normalizeLikelihoodByEventNumber_
Definition: MuScleFitUtils.h:240
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
histoStyle.weight
int weight
Definition: histoStyle.py:51
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 alignCSCRings::e, siStripFEDMonitor_P5_cff::Max, and Pi.

Referenced by MuScleFitUtils::fitMass().

56  {
57  return (0.5 * par[0] * par[1] / TMath::Pi()) /
58  TMath::Max(1.e-10, (x[0] - par[2]) * (x[0] - par[2]) + .25 * par[1] * par[1]);
59 }
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
siStripFEDMonitor_P5_cff.Max
tuple Max
Definition: siStripFEDMonitor_P5_cff.py:56
gpuClustering::x
uint16_t const *__restrict__ x
Definition: gpuClustering.h:43

Variable Documentation

double f[11][100]

Definition at line 78 of file MuScleFitUtils.cc.

double g[11][100]

Definition at line 79 of file MuScleFitUtils.cc.

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().

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.

double isum

Definition at line 77 of file MuScleFitUtils.cc.

Referenced by OuterTrackerMonitorTrackingParticles::analyze(), L1TrackObjectNtupleMaker::analyze(), l1t::CaloParamsHelperO2O::etSumEtaMax(), l1t::CaloParamsHelper::etSumEtaMax(), l1t::CaloParamsHelperO2O::etSumEtaMin(), l1t::CaloParamsHelper::etSumEtaMin(), l1t::CaloParamsHelperO2O::etSumEtThreshold(), l1t::CaloParamsHelper::etSumEtThreshold(), reco::HcalNoiseInfoProducer::filldigis(), SiPixelActionExecutor::fillFEDErrorSummary(), SiStripSummaryCreator::fillGrandSummaryHistos(), SiPixelActionExecutor::fillSummary(), SiStripSummaryCreator::fillSummaryHistos(), main(), CastorDigiMonitor::processEvent(), l1t::CaloParamsHelperO2O::setEtSumEtaMax(), l1t::CaloParamsHelper::setEtSumEtaMax(), l1t::CaloParamsHelperO2O::setEtSumEtaMin(), l1t::CaloParamsHelper::setEtSumEtaMin(), l1t::CaloParamsHelperO2O::setEtSumEtThreshold(), l1t::CaloParamsHelper::setEtSumEtThreshold(), and SiStripSummaryCreator::setSummaryMENames().

double mzsum

Definition at line 76 of file MuScleFitUtils.cc.

Generated for CMSSW Reference Manual by doxygen 1.8.5