#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

Double_t Gaussian	(	Double_t *	x,
		Double_t *	par
	)

Definition at line 63 of file MuScleFitUtils.cc.

References JetChargeProducer_cfi::exp.

Referenced by MuScleFitUtils::fitReso().

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

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

Definition at line 1709 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_, cmsBatch::log, MuScleFitUtils::loopCounter, ResonanceBuilder::mass, MuScleFitUtils::massProb(), MuScleFitUtils::massResolution(), MuScleFitUtils::minuitLoop_, MuScleFitUtils::normalizationChanged_, MuScleFitUtils::normalizeLikelihoodByEventNumber_, MuScleFitUtils::oldNormalization_, TtFullHadEvtBuilder_cfi::prob, MuScleFitUtils::ReducedSavedPair, MuScleFitUtils::rminPtr_, MuScleFitUtils::SavedPair, mps_merge::weight, and DOFs::Y.

Referenced by MuScleFitUtils::minimizeLikelihood(), MuScleFitUtils::byPt::operator()(), 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
 }

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, Max(), and Pi.

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

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 87 of file MuScleFitUtils.cc.

double isum

Definition at line 77 of file MuScleFitUtils.cc.

double mzsum

Definition at line 76 of file MuScleFitUtils.cc.

Functions

Variables

Function Documentation

Variable Documentation