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

typedef reco::Particle::LorentzVector lorentzVector

Definition at line 47 of file MuScleFitUtils.h.

Function Documentation

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

                                                                                          {
  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
}