|
|
|
#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 <vector>Go to the source code of this file.
Classes | |
| struct | MuScleFitUtils::byPt |
| struct | MuScleFitUtils::massResolComponentsStruct |
| class | MuScleFitUtils |
Typedefs | |
| typedef reco::Particle::LorentzVector | lorentzVector |
Functions | |
| void | likelihood (int &npar, double *grad, double &fval, double *xval, int flag) |
Definition at line 42 of file MuScleFitUtils.h.
| void likelihood | ( | int & | npar, |
| double * | grad, | ||
| double & | fval, | ||
| double * | xval, | ||
| int | flag | ||
| ) |
Definition at line 1648 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_, funct::log(), MuScleFitUtils::loopCounter, 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().
{
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 );
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_ == 0 ) {
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_ != 0 ) {
++MuScleFitUtils::minuitLoop_;
MuScleFitUtils::likelihoodInLoop_->SetBinContent(MuScleFitUtils::minuitLoop_, fval);
}
}
// else std::cout << "minuitLoop over 10000. Not filling histogram" << 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
}
1.7.3