CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Static Public Member Functions | Static Public Attributes
MuScleFitUtils Class Reference

#include <MuScleFitUtils.h>

Classes

struct  byPt
 
struct  massResolComponentsStruct
 

Public Member Functions

 MuScleFitUtils ()
 
virtual ~MuScleFitUtils ()
 

Static Public Member Functions

static lorentzVector applyBias (const lorentzVector &muon, const int charge)
 
static lorentzVector applyScale (const lorentzVector &muon, const std::vector< double > &parval, const int charge)
 
static lorentzVector applyScale (const lorentzVector &muon, double *parval, const int charge)
 
static lorentzVector applySmearing (const lorentzVector &muon)
 
static bool checkMassWindow (const double &mass, const double &leftBorder, const double &rightBorder)
 Method to check if the mass value is within the mass window of the i-th resonance. More...
 
static double computeWeight (const double &mass, const int iev, const bool doUseBkgrWindow=false)
 
static double deltaPhi (const double &phi1, const double &phi2)
 
static double deltaPhiNoFabs (const double &phi1, const double &phi2)
 Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance computations. More...
 
static double deltaR (const double &eta1, const double &eta2, const double &phi1, const double &phi2)
 
static std::pair< MuScleFitMuon, MuScleFitMuonfindBestRecoRes (const std::vector< MuScleFitMuon > &muons)
 
static std::pair< SimTrack, SimTrackfindBestSimuRes (const std::vector< SimTrack > &simMuons)
 
static std::pair< lorentzVector, lorentzVectorfindGenMuFromRes (const reco::GenParticleCollection *genParticles)
 
static std::pair< lorentzVector, lorentzVectorfindGenMuFromRes (const edm::HepMCProduct *evtMC)
 
static std::pair< lorentzVector, lorentzVectorfindSimMuFromRes (const edm::Handle< edm::HepMCProduct > &evtMC, const edm::Handle< edm::SimTrackContainer > &simTracks)
 
static std::vector< TGraphErrors * > fitMass (TH2F *histo)
 
static std::vector< TGraphErrors * > fitReso (TH2F *histo)
 
static lorentzVector fromPtEtaPhiToPxPyPz (const double *ptEtaPhiE)
 
static double invDimuonMass (const lorentzVector &mu1, const lorentzVector &mu2)
 
static double massProb (const double &mass, const double &rapidity, const int ires, const double &massResol)
 
static double massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, const std::vector< double > &parval, const bool doUseBkgrWindow, const double &eta1, const double &eta2)
 
static double massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, double *parval, const bool doUseBkgrWindow, const double &eta1, const double &eta2)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, std::unique_ptr< double > parval)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, double *parval)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, const ResolutionFunction &resolFunc)
 
static void minimizeLikelihood ()
 
static double probability (const double &mass, const double &massResol, const double GLvalue[][1001][1001], const double GLnorm[][1001], const int iRes, const int iY)
 Computes the probability given the mass, mass resolution and the arrays with the probabilities and the normalizations. More...
 

Static Public Attributes

static const int backgroundFunctionsRegions
 
static BackgroundHandlerbackgroundHandler
 
static TH1D * backgroundProb_ = nullptr
 
static int BgrFitType = 0
 
static scaleFunctionBase< std::vector< double > > * biasFunction = nullptr
 
static int BiasType = 0
 
static bool computeMinosErrors_
 
static int counter_resprob = 0
 
static double crossSection [6]
 
static CrossSectionHandlercrossSectionHandler
 
static int debug = 0
 
static bool debugMassResol_
 
static double deltaPhiMaxCut_ = 100.
 
static double deltaPhiMinCut_ = -100.
 
static std::vector< int > doBackgroundFit
 
static std::vector< int > doCrossSectionFit
 
static std::vector< int > doResolFit
 
static std::vector< int > doScaleFit
 
static bool duringMinos_ = false
 
static int FitStrategy = 1
 
static std::vector< std::pair< MuScleFitMuon, MuScleFitMuon > > genMuscleFitPair
 
static std::vector< std::pair< lorentzVector, lorentzVector > > genPair
 
static double GLNorm [6][1001]
 
static double GLValue [6][1001][1001]
 
static double GLZNorm [40][1001]
 
static double GLZValue [40][1001][1001]
 
static int goodmuon = 0
 
static int iev_ = 0
 
static TH1D * likelihoodInLoop_ = nullptr
 
static unsigned int loopCounter = 5
 
static struct MuScleFitUtils::massResolComponentsStruct massResolComponents
 
static double massWindowHalfWidth [3][6]
 
static double maxMuonEtaFirstRange_ = 6.
 
static double maxMuonEtaSecondRange_ = 100.
 
static double maxMuonPt_ = 100000000.
 
static bool minimumShapePlots_
 
static double minMuonEtaFirstRange_ = -6.
 
static double minMuonEtaSecondRange_ = -100.
 
static double minMuonPt_ = 0.
 
static int minuitLoop_ = 0
 
static const double mMu2 = 0.011163612
 
static const unsigned int motherPdgIdArray [6] = {23, 100553, 100553, 553, 100443, 443}
 
static const double muMass = 0.105658
 
static int MuonType
 
static int MuonTypeForCheckMassWindow
 
static int nbins = 1000
 
static unsigned int normalizationChanged_ = 0
 
static bool normalizeLikelihoodByEventNumber_ = true
 
static double oldNormalization_ = 0.
 
static std::vector< double > parBgr
 
static std::vector< int > parBgrFix
 
static std::vector< int > parBgrOrder
 
static std::vector< double > parBias
 
static std::vector< double > parCrossSection
 
static std::vector< int > parCrossSectionFix
 
static std::vector< int > parCrossSectionOrder
 
static std::vector< int > parfix
 
static std::vector< int > parorder
 
static std::vector< double > parResol
 
static std::vector< int > parResolFix
 
static std::vector< double > parResolMax
 
static std::vector< double > parResolMin
 
static std::vector< int > parResolOrder
 
static std::vector< double > parResolStep
 
static std::vector< double > parScale
 
static std::vector< int > parScaleFix
 
static std::vector< double > parScaleMax
 
static std::vector< double > parScaleMin
 
static std::vector< int > parScaleOrder
 
static std::vector< double > parScaleStep
 
static std::vector< double > parSmear
 
static std::vector< std::vector< double > > parvalue
 
static bool rapidityBinsForZ_ = true
 
static std::vector< std::pair< lorentzVector, lorentzVector > > ReducedSavedPair
 
static std::vector< int > resfind
 
static bool ResFound = false
 
static double ResGamma [6] = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932}
 
static double ResHalfWidth [6]
 
static double ResMass [6] = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969}
 
static double ResMaxSigma [6]
 
static double ResMinMass [6] = {-99, -99, -99, -99, -99, -99}
 
static int ResolFitType = 0
 
static resolutionFunctionBase< double * > * resolutionFunction = nullptr
 
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec = nullptr
 
static TMinuit * rminPtr_ = nullptr
 
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
 
static std::vector< std::pair< MuScleFitMuon, MuScleFitMuon > > SavedPairMuScleFitMuons
 
static bool scaleFitNotDone_ = true
 
static int ScaleFitType = 0
 
static scaleFunctionBase< double * > * scaleFunction = nullptr
 
static scaleFunctionBase< std::vector< double > > * scaleFunctionForVec = nullptr
 
static bool separateRanges_ = true
 
static bool sherpa_ = false
 
static TH1D * signalProb_ = nullptr
 
static std::vector< std::pair< lorentzVector, lorentzVector > > simPair
 
static smearFunctionBasesmearFunction = nullptr
 
static int SmearType = 0
 
static bool speedup = false
 
static bool startWithSimplex_
 
static const int totalResNum = 6
 
static bool useProbsFile_ = true
 
static double x [7][10000]
 

Detailed Description

Definition at line 51 of file MuScleFitUtils.h.

Constructor & Destructor Documentation

◆ MuScleFitUtils()

MuScleFitUtils::MuScleFitUtils ( )
inline

Definition at line 55 of file MuScleFitUtils.h.

55 {}

◆ ~MuScleFitUtils()

virtual MuScleFitUtils::~MuScleFitUtils ( )
inlinevirtual

Definition at line 59 of file MuScleFitUtils.h.

59 {}

Member Function Documentation

◆ applyBias()

lorentzVector MuScleFitUtils::applyBias ( const lorentzVector muon,
const int  charge 
)
static

Definition at line 444 of file MuScleFitUtils.cc.

References biasFunction, chg, gather_cfg::cout, debug, fromPtEtaPhiToPxPyPz(), parBias, and scaleFunctionBase< T >::scale().

Referenced by MuScleFit::applyBias().

444  {
445  double ptEtaPhiE[4] = {muon.Pt(), muon.Eta(), muon.Phi(), muon.E()};
446 
447  if (MuScleFitUtils::debug > 1)
448  std::cout << "pt before bias = " << ptEtaPhiE[0] << std::endl;
449 
450  // Use functors (although not with the () operator)
451  // Note that we always pass pt, eta and phi, but internally only the needed
452  // values are used.
453  // The functors used are takend from the same group used for the scaling
454  // thus the name of the method used is "scale".
455  ptEtaPhiE[0] = biasFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, MuScleFitUtils::parBias);
456 
457  if (MuScleFitUtils::debug > 1)
458  std::cout << "pt after bias = " << ptEtaPhiE[0] << std::endl;
459 
460  return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
461 }
static std::vector< double > parBias
const float chg[109]
Definition: CoreSimTrack.cc:5
static int debug
static scaleFunctionBase< std::vector< double > > * biasFunction
static lorentzVector fromPtEtaPhiToPxPyPz(const double *ptEtaPhiE)
virtual double scale(const double &pt, const double &eta, const double &phi, const int chg, const T &parScale) const =0

◆ applyScale() [1/2]

lorentzVector MuScleFitUtils::applyScale ( const lorentzVector muon,
const std::vector< double > &  parval,
const int  charge 
)
static

Definition at line 465 of file MuScleFitUtils.cc.

References chg, EcalPhiSymFlatTableProducers_cfi::id, createfilelist::int, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and AlCaHLTBitMon_ParallelJobs::p.

Referenced by MuScleFit::duringFastLoop(), and likelihood().

465  {
466  double *p = new double[(int)(parval.size())];
467  // Replaced by auto_ptr, which handles delete at the end
468  // std::unique_ptr<double> p(new double[(int)(parval.size())]);
469  // Removed auto_ptr, check massResolution for an explanation.
470  int id = 0;
471  for (std::vector<double>::const_iterator it = parval.begin(); it != parval.end(); ++it, ++id) {
472  //(&*p)[id] = *it;
473  // Also ok would be (p.get())[id] = *it;
474  p[id] = *it;
475  }
476  lorentzVector tempScaleVec(applyScale(muon, p, chg));
477  delete[] p;
478  return tempScaleVec;
479 }
const float chg[109]
Definition: CoreSimTrack.cc:5
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)

◆ applyScale() [2/2]

lorentzVector MuScleFitUtils::applyScale ( const lorentzVector muon,
double *  parval,
const int  charge 
)
static

Definition at line 483 of file MuScleFitUtils.cc.

References chg, gather_cfg::cout, debug, fromPtEtaPhiToPxPyPz(), parResol, scaleFunctionBase< T >::scale(), scaleFunction, and edm::shift.

483  {
484  double ptEtaPhiE[4] = {muon.Pt(), muon.Eta(), muon.Phi(), muon.E()};
485  int shift = parResol.size();
486 
487  if (MuScleFitUtils::debug > 1)
488  std::cout << "pt before scale = " << ptEtaPhiE[0] << std::endl;
489 
490  // the address of parval[shift] is passed as pointer to double. Internally it is used as a normal array, thus:
491  // array[0] = parval[shift], array[1] = parval[shift+1], ...
492  ptEtaPhiE[0] = scaleFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, &(parval[shift]));
493 
494  if (MuScleFitUtils::debug > 1)
495  std::cout << "pt after scale = " << ptEtaPhiE[0] << std::endl;
496 
497  return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
498 }
const float chg[109]
Definition: CoreSimTrack.cc:5
static std::vector< double > parResol
static int debug
static lorentzVector fromPtEtaPhiToPxPyPz(const double *ptEtaPhiE)
virtual double scale(const double &pt, const double &eta, const double &phi, const int chg, const T &parScale) const =0
static unsigned int const shift
static scaleFunctionBase< double * > * scaleFunction

◆ applySmearing()

lorentzVector MuScleFitUtils::applySmearing ( const lorentzVector muon)
static

Definition at line 416 of file MuScleFitUtils.cc.

References gather_cfg::cout, debug, PVValHelper::eta, fromPtEtaPhiToPxPyPz(), goodmuon, mps_fire::i, parSmear, phi, DiDispStaMuonMonitor_cfi::pt, smearFunctionBase::smear(), smearFunction, SmearType, x, and y.

Referenced by MuScleFit::applySmearing().

416  {
417  double pt = muon.Pt();
418  double eta = muon.Eta();
419  double phi = muon.Phi();
420  double E = muon.E();
421 
422  double y[7] = {};
423  for (int i = 0; i < SmearType + 3; i++) {
424  y[i] = x[i][goodmuon % 10000];
425  }
426 
427  // Use the smear function selected in the constructor
429 
430  if (debug > 9) {
431  std::cout << "Smearing Pt,eta,phi = " << pt << " " << eta << " " << phi << "; x = ";
432  for (int i = 0; i < SmearType + 3; i++) {
433  std::cout << y[i];
434  }
435  std::cout << std::endl;
436  }
437 
438  double ptEtaPhiE[4] = {pt, eta, phi, E};
439  return (fromPtEtaPhiToPxPyPz(ptEtaPhiE));
440 }
static smearFunctionBase * smearFunction
virtual void smear(double &pt, double &eta, double &phi, const double *y, const std::vector< double > &parSmear)=0
static double x[7][10000]
static int debug
static int SmearType
static lorentzVector fromPtEtaPhiToPxPyPz(const double *ptEtaPhiE)
static std::vector< double > parSmear
static int goodmuon

◆ checkMassWindow()

bool MuScleFitUtils::checkMassWindow ( const double &  mass,
const double &  leftBorder,
const double &  rightBorder 
)
inlinestatic

Method to check if the mass value is within the mass window of the i-th resonance.

Definition at line 1157 of file MuScleFitUtils.cc.

References EgHLTOffHistBins_cfi::mass.

Referenced by computeWeight(), massProb(), and minimizeLikelihood().

1157  {
1158  return ((mass > leftBorder) && (mass < rightBorder));
1159 }

◆ computeWeight()

double MuScleFitUtils::computeWeight ( const double &  mass,
const int  iev,
const bool  doUseBkgrWindow = false 
)
static

Definition at line 1163 of file MuScleFitUtils.cc.

References backgroundHandler, checkMassWindow(), gather_cfg::cout, debug, doBackgroundFit, loopCounter, EgHLTOffHistBins_cfi::mass, resfind, mps_merge::weight, and BackgroundHandler::windowBorders().

Referenced by MuScleFit::duringFastLoop(), and likelihood().

1163  {
1164  // Compute weight for this event
1165  // -----------------------------
1166  double weight = 0.;
1167 
1168  // Take the highest-mass resonance within bounds
1169  // NB this must be revised once credible estimates of the relative xs of Y(1S), (2S), and (3S)
1170  // are made. Those are priors in the decision of which resonance to assign to an in-between event.
1171  // -----------------------------------------------------------------------------------------------
1172 
1173  if (doUseBkgrWindow && (debug > 0))
1174  std::cout << "using backgrond window for mass = " << mass << std::endl;
1175  // bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
1176  bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
1177 
1178  for (int ires = 0; ires < 6; ires++) {
1179  if (resfind[ires] > 0 && weight == 0.) {
1180  // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
1181  std::pair<double, double> windowBorder = backgroundHandler->windowBorders(useBackgroundWindow, ires);
1182  // if( checkMassWindow(mass, ires, backgroundHandler->resMass( useBackgroundWindow, ires ),
1183  // windowFactor.first, windowFactor.second) ) {
1184  if (checkMassWindow(mass, windowBorder.first, windowBorder.second)) {
1185  weight = 1.0;
1186  if (doUseBkgrWindow && (debug > 0))
1187  std::cout << "setting weight to = " << weight << std::endl;
1188  }
1189  }
1190  }
1191 
1192  return weight;
1193 }
static unsigned int loopCounter
static std::vector< int > doBackgroundFit
static BackgroundHandler * backgroundHandler
static int debug
Definition: weight.py:1
static bool checkMassWindow(const double &mass, const double &leftBorder, const double &rightBorder)
Method to check if the mass value is within the mass window of the i-th resonance.
std::pair< double, double > windowBorders(const bool doBackgroundFit, const int ires)
Returns the appropriate window borders depending on whether the background is being fitted and on the...
static std::vector< int > resfind

◆ deltaPhi()

static double MuScleFitUtils::deltaPhi ( const double &  phi1,
const double &  phi2 
)
inlinestatic

Definition at line 109 of file MuScleFitUtils.h.

References Pi.

Referenced by ResolutionAnalyzer::checkDeltaR(), MuScleFit::checkDeltaR(), deltaPhiNoFabs(), deltaR(), HDelta::Fill(), HMassResolutionVSPart::Fill(), and MuScleFit::selectMuons().

109  {
110  double deltaPhi = phi1 - phi2;
111  while (deltaPhi >= TMath::Pi())
112  deltaPhi -= 2 * TMath::Pi();
113  while (deltaPhi < -TMath::Pi())
114  deltaPhi += 2 * TMath::Pi();
115  return fabs(deltaPhi);
116  }
const double Pi
static double deltaPhi(const double &phi1, const double &phi2)

◆ deltaPhiNoFabs()

static double MuScleFitUtils::deltaPhiNoFabs ( const double &  phi1,
const double &  phi2 
)
inlinestatic

Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance computations.

Definition at line 118 of file MuScleFitUtils.h.

References deltaPhi(), and Pi.

Referenced by ResolutionAnalyzer::analyze(), HCovarianceVSParts::Fill(), and MuScleFit::fillComparisonHistograms().

118  {
119  double deltaPhi = phi1 - phi2;
120  while (deltaPhi >= TMath::Pi())
121  deltaPhi -= 2 * TMath::Pi();
122  while (deltaPhi < -TMath::Pi())
123  deltaPhi += 2 * TMath::Pi();
124  return deltaPhi;
125  }
const double Pi
static double deltaPhi(const double &phi1, const double &phi2)

◆ deltaR()

static double MuScleFitUtils::deltaR ( const double &  eta1,
const double &  eta2,
const double &  phi1,
const double &  phi2 
)
inlinestatic

Definition at line 126 of file MuScleFitUtils.h.

References deltaPhi(), HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, funct::pow(), and mathSSE::sqrt().

126  {
127  return sqrt(std::pow(eta1 - eta2, 2) + std::pow(deltaPhi(phi1, phi2), 2));
128  }
T sqrt(T t)
Definition: SSEVec.h:23
static double deltaPhi(const double &phi1, const double &phi2)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ findBestRecoRes()

std::pair< MuScleFitMuon, MuScleFitMuon > MuScleFitUtils::findBestRecoRes ( const std::vector< MuScleFitMuon > &  muons)
static

Definition at line 315 of file MuScleFitUtils.cc.

References funct::abs(), gather_cfg::cout, debug, HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, EgHLTOffHistBins_cfi::mass, massProb(), massResolution(), maxMuonEtaFirstRange_, maxMuonPt_, minMuonEtaFirstRange_, minMuonEtaSecondRange_, minMuonPt_, HLT_2024v14_cff::Muon1, DiMuonV_cfg::muons, parResol, TtFullHadEvtBuilder_cfi::prob, HLT_2024v14_cff::pt1, HLT_2024v14_cff::pt2, resfind, ResFound, ResMass, useProbsFile_, and beamSpotPI::Y.

Referenced by TestCorrection::analyze(), and MuScleFit::selectMuons().

315  {
316  // NB this routine returns the resonance, but it also sets the ResFound flag, which
317  // is used in MuScleFit to decide whether to use the event or not.
318  // --------------------------------------------------------------------------------
319  if (debug > 0)
320  std::cout << "In findBestRecoRes" << std::endl;
321  ResFound = false;
322  std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes;
323 
324  // Choose the best resonance using its mass probability
325  // ----------------------------------------------------
326  double maxprob = -0.1;
327  double minDeltaMass = 999999;
328  std::pair<MuScleFitMuon, MuScleFitMuon> bestMassMuons;
329  for (std::vector<MuScleFitMuon>::const_iterator Muon1 = muons.begin(); Muon1 != muons.end(); ++Muon1) {
330  //rc2010
331  if (debug > 0)
332  std::cout << "muon_1_charge:" << (*Muon1).charge() << std::endl;
333  for (std::vector<MuScleFitMuon>::const_iterator Muon2 = Muon1 + 1; Muon2 != muons.end(); ++Muon2) {
334  //rc2010
335  if (debug > 0)
336  std::cout << "after_2" << std::endl;
337  if ((((*Muon1).charge()) * ((*Muon2).charge())) > 0) {
338  continue; // This also gets rid of Muon1==Muon2...
339  }
340  // To allow the selection of ranges at negative and positive eta independently we define two
341  // ranges of eta: (minMuonEtaFirstRange_, maxMuonEtaFirstRange_) and (minMuonEtaSecondRange_, maxMuonEtaSecondRange_).
342  double ch1 = (*Muon1).charge();
343  double ch2 = (*Muon2).charge();
344  double pt1 = (*Muon1).Pt();
345  double pt2 = (*Muon2).Pt();
346  double eta1 = (*Muon1).Eta();
347  double eta2 = (*Muon2).Eta();
348  if (pt1 >= minMuonPt_ && pt1 < maxMuonPt_ && pt2 >= minMuonPt_ && pt2 < maxMuonPt_ &&
349  ((eta1 >= minMuonEtaFirstRange_ && eta1 < maxMuonEtaFirstRange_ && eta2 >= minMuonEtaSecondRange_ &&
350  eta2 < maxMuonEtaSecondRange_ && ch1 >= ch2 // In the configuration file, MuonOne==MuonPlus
351  ) ||
352  (eta1 >= minMuonEtaSecondRange_ && eta1 < maxMuonEtaSecondRange_ && eta2 >= minMuonEtaFirstRange_ &&
353  eta2 < maxMuonEtaFirstRange_ && ch1 < ch2))) {
354  double mcomb = ((*Muon1).p4() + (*Muon2).p4()).mass();
355  double Y = ((*Muon1).p4() + (*Muon2).p4()).Rapidity();
356  if (debug > 1) {
357  std::cout << "muon1 " << (*Muon1).p4().Px() << ", " << (*Muon1).p4().Py() << ", " << (*Muon1).p4().Pz()
358  << ", " << (*Muon1).p4().E() << ", " << (*Muon1).charge() << std::endl;
359  std::cout << "muon2 " << (*Muon2).p4().Px() << ", " << (*Muon2).p4().Py() << ", " << (*Muon2).p4().Pz()
360  << ", " << (*Muon2).p4().E() << ", " << (*Muon2).charge() << std::endl;
361  std::cout << "mcomb " << mcomb << std::endl;
362  }
363  double massResol = 0.;
364  if (useProbsFile_) {
365  massResol = massResolution((*Muon1).p4(), (*Muon2).p4(), parResol);
366  }
367  double prob = 0;
368  for (int ires = 0; ires < 6; ires++) {
369  if (resfind[ires] > 0) {
370  if (useProbsFile_) {
371  prob = massProb(mcomb, Y, ires, massResol);
372  }
373  if (prob > maxprob) {
374  if (ch1 < 0) { // store first the mu minus and then the mu plus
375  recMuFromBestRes.first = (*Muon1);
376  recMuFromBestRes.second = (*Muon2);
377  } else {
378  recMuFromBestRes.first = (*Muon2);
379  recMuFromBestRes.second = (*Muon1);
380  }
381  if (debug > 0)
382  std::cout << "muon_1_charge (after swapping):" << recMuFromBestRes.first.charge() << std::endl;
383  ResFound = true; // NNBB we accept "resonances" even outside mass bounds
384  maxprob = prob;
385  }
386  // if( ResMass[ires] == 0 ) {
387  // std::cout << "Error: ResMass["<<ires<<"] = " << ResMass[ires] << std::endl;
388  // exit(1);
389  // }
390  double deltaMass = std::abs(mcomb - ResMass[ires]) / ResMass[ires];
391  if (deltaMass < minDeltaMass) {
392  bestMassMuons = std::make_pair((*Muon1), (*Muon2));
393  minDeltaMass = deltaMass;
394  }
395  }
396  }
397  }
398  }
399  }
400  //If outside mass window (maxprob==0) then take the two muons with best invariant mass
401  //(anyway they will not be used in the likelihood calculation, only to fill plots)
402  if (!maxprob) {
403  if (bestMassMuons.first.charge() < 0) {
404  recMuFromBestRes.first = bestMassMuons.first;
405  recMuFromBestRes.second = bestMassMuons.second;
406  } else {
407  recMuFromBestRes.second = bestMassMuons.first;
408  recMuFromBestRes.first = bestMassMuons.second;
409  }
410  }
411  return recMuFromBestRes;
412 }
static std::vector< double > parResol
static int debug
static double ResMass[6]
static bool ResFound
static double minMuonEtaFirstRange_
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
muons
the two sets of parameters below are mutually exclusive, depending if RECO or ALCARECO is used the us...
Definition: DiMuonV_cfg.py:214
static double maxMuonPt_
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
static double minMuonEtaSecondRange_
static double minMuonPt_
static bool useProbsFile_
static std::vector< int > resfind
static double maxMuonEtaFirstRange_

◆ findBestSimuRes()

std::pair< SimTrack, SimTrack > MuScleFitUtils::findBestSimuRes ( const std::vector< SimTrack > &  simMuons)
static

Definition at line 279 of file MuScleFitUtils.cc.

References EgHLTOffHistBins_cfi::mass, massProb(), TtFullHadEvtBuilder_cfi::prob, resfind, and beamSpotPI::Y.

Referenced by MuScleFitPlotter::fillSim().

279  {
280  std::pair<SimTrack, SimTrack> simMuFromBestRes;
281  double maxprob = -0.1;
282 
283  // Double loop on muons
284  // --------------------
285  for (std::vector<SimTrack>::const_iterator simMu1 = simMuons.begin(); simMu1 != simMuons.end(); simMu1++) {
286  for (std::vector<SimTrack>::const_iterator simMu2 = simMu1 + 1; simMu2 != simMuons.end(); simMu2++) {
287  if (((*simMu1).charge() * (*simMu2).charge()) > 0) {
288  continue; // this also gets rid of simMu1==simMu2...
289  }
290  // Choose the best resonance using its mass. Check Z, Y(3S,2S,1S), Psi(2S), J/Psi in order
291  // ---------------------------------------------------------------------------------------
292  double mcomb = ((*simMu1).momentum() + (*simMu2).momentum()).mass();
293  double Y = ((*simMu1).momentum() + (*simMu2).momentum()).Rapidity();
294  for (int ires = 0; ires < 6; ires++) {
295  if (resfind[ires] > 0) {
296  double prob = massProb(mcomb, Y, ires, 0.);
297  if (prob > maxprob) {
298  simMuFromBestRes.first = (*simMu1);
299  simMuFromBestRes.second = (*simMu2);
300  maxprob = prob;
301  }
302  }
303  }
304  }
305  }
306 
307  // Return most likely combination of muons making a resonance
308  // ----------------------------------------------------------
309  return simMuFromBestRes;
310 }
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
static std::vector< int > resfind

◆ findGenMuFromRes() [1/2]

std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findGenMuFromRes ( const reco::GenParticleCollection genParticles)
static

Definition at line 2344 of file MuScleFitUtils.cc.

References funct::abs(), gather_cfg::cout, debug, AJJGenJetFilter_cfi::genParticles, motherPdgIdArray, and resfind.

Referenced by MuScleFitGenFilter::filter().

2345  {
2346  std::pair<lorentzVector, lorentzVector> muFromRes;
2347 
2348  //Loop on generated particles
2349  if (debug > 0)
2350  std::cout << "Starting loop on " << genParticles->size() << " genParticles" << std::endl;
2351  for (reco::GenParticleCollection::const_iterator part = genParticles->begin(); part != genParticles->end(); ++part) {
2352  if (std::abs(part->pdgId()) == 13 && part->status() == 1) {
2353  bool fromRes = false;
2354  unsigned int motherPdgId = part->mother()->pdgId();
2355  if (debug > 0) {
2356  std::cout << "Found a muon with mother: " << motherPdgId << std::endl;
2357  }
2358  for (int ires = 0; ires < 6; ++ires) {
2359  if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2360  fromRes = true;
2361  }
2362  if (fromRes) {
2363  if (part->pdgId() == 13) {
2364  muFromRes.first = part->p4();
2365  if (debug > 0)
2366  std::cout << "Found a genMuon + : " << muFromRes.first << std::endl;
2367  // muFromRes.first = (lorentzVector(part->p4().px(),part->p4().py(),
2368  // part->p4().pz(),part->p4().e()));
2369  } else {
2370  muFromRes.second = part->p4();
2371  if (debug > 0)
2372  std::cout << "Found a genMuon - : " << muFromRes.second << std::endl;
2373  // muFromRes.second = (lorentzVector(part->p4().px(),part->p4().py(),
2374  // part->p4().pz(),part->p4().e()));
2375  }
2376  }
2377  }
2378  }
2379  return muFromRes;
2380 }
static int debug
static const unsigned int motherPdgIdArray[6]
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
part
Definition: HCALResponse.h:20
static std::vector< int > resfind

◆ findGenMuFromRes() [2/2]

std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findGenMuFromRes ( const edm::HepMCProduct evtMC)
static

Definition at line 2306 of file MuScleFitUtils.cc.

References funct::abs(), edm::HepMCProduct::GetEvent(), motherPdgIdArray, resfind, and sherpa_.

2306  {
2307  const HepMC::GenEvent *Evt = evtMC->GetEvent();
2308  std::pair<lorentzVector, lorentzVector> muFromRes;
2309  //Loop on generated particles
2310  for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); part++) {
2311  if (std::abs((*part)->pdg_id()) == 13 && (*part)->status() == 1) {
2312  bool fromRes = false;
2313  for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);
2314  mother != (*part)->production_vertex()->particles_end(HepMC::ancestors);
2315  ++mother) {
2316  unsigned int motherPdgId = (*mother)->pdg_id();
2317 
2318  // For sherpa the resonance is not saved. The muons from the resonance can be identified
2319  // by having as mother a muon of status 3.
2320  if (sherpa_) {
2321  if (motherPdgId == 13 && (*mother)->status() == 3)
2322  fromRes = true;
2323  } else {
2324  for (int ires = 0; ires < 6; ++ires) {
2325  if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2326  fromRes = true;
2327  }
2328  }
2329  }
2330  if (fromRes) {
2331  if ((*part)->pdg_id() == 13)
2332  // muFromRes.first = (*part)->momentum();
2333  muFromRes.first = (lorentzVector(
2334  (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
2335  else
2336  muFromRes.second = (lorentzVector(
2337  (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
2338  }
2339  }
2340  }
2341  return muFromRes;
2342 }
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static const unsigned int motherPdgIdArray[6]
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
static bool sherpa_
part
Definition: HCALResponse.h:20
static std::vector< int > resfind

◆ findSimMuFromRes()

std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findSimMuFromRes ( const edm::Handle< edm::HepMCProduct > &  evtMC,
const edm::Handle< edm::SimTrackContainer > &  simTracks 
)
static

Definition at line 2265 of file MuScleFitUtils.cc.

References funct::abs(), GenParticle::GenParticle, edm::HepMCProduct::GetEvent(), runTauDisplay::gp, motherPdgIdArray, resfind, cscDigiValidation_cfi::simTrack, and TrackCandidateProducer_cfi::simTracks.

Referenced by MuScleFitPlotter::fillGenSim().

2266  {
2267  //Loop on simulated tracks
2268  std::pair<lorentzVector, lorentzVector> simMuFromRes;
2269  for (edm::SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
2270  //Chose muons
2271  if (std::abs((*simTrack).type()) == 13) {
2272  //If tracks from IP than find mother
2273  if ((*simTrack).genpartIndex() > 0) {
2274  HepMC::GenParticle *gp = evtMC->GetEvent()->barcode_to_particle((*simTrack).genpartIndex());
2275  if (gp != nullptr) {
2276  for (HepMC::GenVertex::particle_iterator mother = gp->production_vertex()->particles_begin(HepMC::ancestors);
2277  mother != gp->production_vertex()->particles_end(HepMC::ancestors);
2278  ++mother) {
2279  bool fromRes = false;
2280  unsigned int motherPdgId = (*mother)->pdg_id();
2281  for (int ires = 0; ires < 6; ++ires) {
2282  if (motherPdgId == motherPdgIdArray[ires] && resfind[ires])
2283  fromRes = true;
2284  }
2285  if (fromRes) {
2286  if (gp->pdg_id() == 13)
2287  simMuFromRes.first = lorentzVector(simTrack->momentum().px(),
2288  simTrack->momentum().py(),
2289  simTrack->momentum().pz(),
2290  simTrack->momentum().e());
2291  else
2292  simMuFromRes.second = lorentzVector(simTrack->momentum().px(),
2293  simTrack->momentum().py(),
2294  simTrack->momentum().pz(),
2295  simTrack->momentum().e());
2296  }
2297  }
2298  }
2299  // else LogDebug("MuScleFitUtils") << "WARNING: no matching genParticle found for simTrack" << std::endl;
2300  }
2301  }
2302  }
2303  return simMuFromRes;
2304 }
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static const unsigned int motherPdgIdArray[6]
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
static std::vector< int > resfind

◆ fitMass()

std::vector< TGraphErrors * > MuScleFitUtils::fitMass ( TH2F *  histo)
static

Definition at line 1985 of file MuScleFitUtils.cc.

References isoTrack_cff::chi2, generateEDF::cont, gather_cfg::cout, debug, MillePedeFileConverter_cfg::e, submitPVResolutionJobs::err, timingPdfMaker::histo, mps_fire::i, dqmiolumiharvest::j, lorentzianPeak(), mergeVDriftHistosByStation::name, groupFilesInBlocks::nn, mysort::results, x, and geometryCSVtoXML::xx.

1985  {
1986  if (MuScleFitUtils::debug > 0)
1987  std::cout << "Fitting " << histo->GetName() << std::endl;
1988 
1989  std::vector<TGraphErrors *> results;
1990 
1991  // Results of the fit
1992  // ------------------
1993  std::vector<double> Ftop;
1994  std::vector<double> Fwidth;
1995  std::vector<double> Fmass;
1996  std::vector<double> Etop;
1997  std::vector<double> Ewidth;
1998  std::vector<double> Emass;
1999  std::vector<double> Fchi2;
2000  // X bin center and width
2001  // ----------------------
2002  std::vector<double> Xcenter;
2003  std::vector<double> Ex;
2004 
2005  // Fit with lorentzian peak
2006  // ------------------------
2007  TF1 *fitFcn = new TF1("fitFcn", lorentzianPeak, 70, 110, 3);
2008  fitFcn->SetParameters(100, 3, 91);
2009  fitFcn->SetParNames("Ftop", "Fwidth", "Fmass");
2010  fitFcn->SetLineWidth(2);
2011 
2012  // Fit slices projected along Y from bins in X
2013  // -------------------------------------------
2014  double cont_min = 20; // Minimum number of entries
2015  Int_t binx = histo->GetXaxis()->GetNbins();
2016  // TFile *f= new TFile("prova.root", "recreate");
2017  // histo->Write();
2018  for (int i = 1; i <= binx; i++) {
2019  TH1 *histoY = histo->ProjectionY("", i, i);
2020  // histoY->Write();
2021  double cont = histoY->GetEntries();
2022  if (cont > cont_min) {
2023  histoY->Fit("fitFcn", "0", "", 70, 110);
2024  double *par = fitFcn->GetParameters();
2025  const double *err = fitFcn->GetParErrors();
2026 
2027  Ftop.push_back(par[0]);
2028  Fwidth.push_back(par[1]);
2029  Fmass.push_back(par[2]);
2030  Etop.push_back(err[0]);
2031  Ewidth.push_back(err[1]);
2032  Emass.push_back(err[2]);
2033 
2034  double chi2 = fitFcn->GetChisquare();
2035  Fchi2.push_back(chi2);
2036 
2037  double xx = histo->GetXaxis()->GetBinCenter(i);
2038  Xcenter.push_back(xx);
2039  double ex = 0; // FIXME: you can use the bin width
2040  Ex.push_back(ex);
2041  }
2042  }
2043  // f->Close();
2044 
2045  // Put the fit results in arrays for TGraphErrors
2046  // ----------------------------------------------
2047  const int nn = Fmass.size();
2048  double *x = new double[nn];
2049  double *ym = new double[nn];
2050  double *e = new double[nn];
2051  double *eym = new double[nn];
2052  double *yw = new double[nn];
2053  double *eyw = new double[nn];
2054  double *yc = new double[nn];
2055 
2056  for (int j = 0; j < nn; j++) {
2057  x[j] = Xcenter[j];
2058  ym[j] = Fmass[j];
2059  eym[j] = Emass[j];
2060  yw[j] = Fwidth[j];
2061  eyw[j] = Ewidth[j];
2062  yc[j] = Fchi2[j];
2063  e[j] = Ex[j];
2064  }
2065 
2066  // Create TGraphErrors
2067  // -------------------
2068  TString name = histo->GetName();
2069  TGraphErrors *grM = new TGraphErrors(nn, x, ym, e, eym);
2070  grM->SetTitle(name + "_M");
2071  grM->SetName(name + "_M");
2072  TGraphErrors *grW = new TGraphErrors(nn, x, yw, e, eyw);
2073  grW->SetTitle(name + "_W");
2074  grW->SetName(name + "_W");
2075  TGraphErrors *grC = new TGraphErrors(nn, x, yc, e, e);
2076  grC->SetTitle(name + "_chi2");
2077  grC->SetName(name + "_chi2");
2078 
2079  // Cleanup
2080  // -------
2081  delete[] x;
2082  delete[] ym;
2083  delete[] eym;
2084  delete[] yw;
2085  delete[] eyw;
2086  delete[] yc;
2087  delete[] e;
2088  delete fitFcn;
2089 
2090  results.push_back(grM);
2091  results.push_back(grW);
2092  results.push_back(grC);
2093  return results;
2094 }
static double x[7][10000]
static int debug
Double_t lorentzianPeak(Double_t *x, Double_t *par)
results
Definition: mysort.py:8
cont
load Luminosity info ##
Definition: generateEDF.py:620

◆ fitReso()

std::vector< TGraphErrors * > MuScleFitUtils::fitReso ( TH2F *  histo)
static

Definition at line 2098 of file MuScleFitUtils.cc.

References isoTrack_cff::chi2, generateEDF::cont, gather_cfg::cout, MillePedeFileConverter_cfg::e, submitPVResolutionJobs::err, Gaussian(), timingPdfMaker::histo, mps_fire::i, dqmiolumiharvest::j, compareTotals::means, mergeVDriftHistosByStation::name, groupFilesInBlocks::nn, mysort::results, x, and geometryCSVtoXML::xx.

2098  {
2099  std::cout << "Fitting " << histo->GetName() << std::endl;
2100  std::vector<TGraphErrors *> results;
2101 
2102  // Results from fit
2103  // ----------------
2104  std::vector<double> maxs;
2105  std::vector<double> means;
2106  std::vector<double> sigmas;
2107  std::vector<double> chi2s;
2108  std::vector<double> Emaxs;
2109  std::vector<double> Emeans;
2110  std::vector<double> Esigmas;
2111 
2112  // X bin center and width
2113  // ----------------------
2114  std::vector<double> Xcenter;
2115  std::vector<double> Ex;
2116 
2117  // Fit with a gaussian
2118  // -------------------
2119  TF1 *fitFcn = new TF1("fitFunc", Gaussian, -0.2, 0.2, 3);
2120  fitFcn->SetParameters(100, 0, 0.02);
2121  fitFcn->SetParNames("max", "mean", "sigma");
2122  fitFcn->SetLineWidth(2);
2123 
2124  // Fit slices projected along Y from bins in X
2125  // -------------------------------------------
2126  double cont_min = 20; // Minimum number of entries
2127  Int_t binx = histo->GetXaxis()->GetNbins();
2128  for (int i = 1; i <= binx; i++) {
2129  TH1 *histoY = histo->ProjectionY("", i, i);
2130  double cont = histoY->GetEntries();
2131  if (cont > cont_min) {
2132  histoY->Fit("fitFunc", "0", "", -0.2, 0.2);
2133  double *par = fitFcn->GetParameters();
2134  const double *err = fitFcn->GetParErrors();
2135 
2136  maxs.push_back(par[0]);
2137  means.push_back(par[1]);
2138  sigmas.push_back(par[2]);
2139  Emaxs.push_back(err[0]);
2140  Emeans.push_back(err[1]);
2141  Esigmas.push_back(err[2]);
2142 
2143  double chi2 = fitFcn->GetChisquare();
2144  chi2s.push_back(chi2);
2145 
2146  double xx = histo->GetXaxis()->GetBinCenter(i);
2147  Xcenter.push_back(xx);
2148  double ex = 0; // FIXME: you can use the bin width
2149  Ex.push_back(ex);
2150  }
2151  }
2152 
2153  // Put the fit results in arrays for TGraphErrors
2154  // ----------------------------------------------
2155  const int nn = means.size();
2156  double *x = new double[nn];
2157  double *ym = new double[nn];
2158  double *e = new double[nn];
2159  double *eym = new double[nn];
2160  double *yw = new double[nn];
2161  double *eyw = new double[nn];
2162  double *yc = new double[nn];
2163 
2164  for (int j = 0; j < nn; j++) {
2165  x[j] = Xcenter[j];
2166  ym[j] = means[j];
2167  eym[j] = Emeans[j];
2168  // yw[j] = maxs[j];
2169  // eyw[j] = Emaxs[j];
2170  yw[j] = sigmas[j];
2171  eyw[j] = Esigmas[j];
2172  yc[j] = chi2s[j];
2173  e[j] = Ex[j];
2174  }
2175 
2176  // Create TGraphErrors
2177  // -------------------
2178  TString name = histo->GetName();
2179  TGraphErrors *grM = new TGraphErrors(nn, x, ym, e, eym);
2180  grM->SetTitle(name + "_mean");
2181  grM->SetName(name + "_mean");
2182  TGraphErrors *grW = new TGraphErrors(nn, x, yw, e, eyw);
2183  grW->SetTitle(name + "_sigma");
2184  grW->SetName(name + "_sigma");
2185  TGraphErrors *grC = new TGraphErrors(nn, x, yc, e, e);
2186  grC->SetTitle(name + "_chi2");
2187  grC->SetName(name + "_chi2");
2188 
2189  // Cleanup
2190  // -------
2191  delete[] x;
2192  delete[] ym;
2193  delete[] eym;
2194  delete[] yw;
2195  delete[] eyw;
2196  delete[] yc;
2197  delete[] e;
2198  delete fitFcn;
2199 
2200  results.push_back(grM);
2201  results.push_back(grW);
2202  results.push_back(grC);
2203  return results;
2204 }
Double_t Gaussian(Double_t *x, Double_t *par)
static double x[7][10000]
results
Definition: mysort.py:8
cont
load Luminosity info ##
Definition: generateEDF.py:620

◆ fromPtEtaPhiToPxPyPz()

lorentzVector MuScleFitUtils::fromPtEtaPhiToPxPyPz ( const double *  ptEtaPhiE)
static

Definition at line 502 of file MuScleFitUtils.cc.

References funct::cos(), JetChargeProducer_cfi::exp, muMass, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, funct::sin(), mathSSE::sqrt(), and createJobs::tmp.

Referenced by applyBias(), applyScale(), applySmearing(), and TestCorrection::correctMuon().

502  {
503  double px = ptEtaPhiE[0] * cos(ptEtaPhiE[2]);
504  double py = ptEtaPhiE[0] * sin(ptEtaPhiE[2]);
505  double tmp = 2 * atan(exp(-ptEtaPhiE[1]));
506  double pz = ptEtaPhiE[0] * cos(tmp) / sin(tmp);
507  double E = sqrt(px * px + py * py + pz * pz + muMass * muMass);
508 
509  // lorentzVector corrMu(px,py,pz,E);
510  // To fix memory leaks, this is to be substituted with
511  // std::unique_ptr<lorentzVector> corrMu(new lorentzVector(px, py, pz, E));
512 
513  return lorentzVector(px, py, pz, E);
514 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static const double muMass
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
tmp
align.sh
Definition: createJobs.py:716

◆ invDimuonMass()

double MuScleFitUtils::invDimuonMass ( const lorentzVector mu1,
const lorentzVector mu2 
)
inlinestatic

Definition at line 518 of file MuScleFitUtils.cc.

Referenced by likelihood(), and minimizeLikelihood().

518  {
519  return (mu1 + mu2).mass();
520 }

◆ massProb() [1/3]

double MuScleFitUtils::massProb ( const double &  mass,
const double &  rapidity,
const int  ires,
const double &  massResol 
)
static

Definition at line 2208 of file MuScleFitUtils.cc.

References gather_cfg::cout, debug, MillePedeFileConverter_cfg::e, GL, EgHLTOffHistBins_cfi::mass, np, Pi, ResGamma, ResMass, w(), and x.

Referenced by MuScleFit::duringFastLoop(), findBestRecoRes(), findBestSimuRes(), likelihood(), and massProb().

2208  {
2209  // This routine computes the likelihood that a given measured mass "measMass" is
2210  // the result of resonance #ires if the resolution expected for the two muons is massResol
2211  // ---------------------------------------------------------------------------------------
2212 
2213  double P = 0.;
2214 
2215  // Return Lorentz value for zero resolution cases (like MC)
2216  // --------------------------------------------------------
2217  if (massResol == 0.) {
2218  if (debug > 9)
2219  std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2220  << massResol << " : used Lorentz P-value" << std::endl;
2221  return (0.5 * ResGamma[ires] / TMath::Pi()) /
2222  ((mass - ResMass[ires]) * (mass - ResMass[ires]) + .25 * ResGamma[ires] * ResGamma[ires]);
2223  }
2224 
2225  // NB defined as below, P is not a "probability" but a likelihood that we observe
2226  // a dimuon mass "mass", given mRef, gamma, and massResol. It is what we need for the
2227  // fit which finds the best resolution parameters, though. A definition which is
2228  // more properly a probability is given below (in massProb2()), where the result
2229  // cannot be used to fit resolution parameters because the fit would always prefer
2230  // to set the res parameters to the minimum possible value (best resolution),
2231  // to have a probability close to one of observing any mass.
2232  // -------------------------------------------------------------------------------
2233  // NNBB the following two lines have been replaced with the six following them,
2234  // which provide an improvement of a factor 9 in speed of execution at a
2235  // negligible price in precision.
2236  // ----------------------------------------------------------------------------
2237  // GL->SetParameters(gamma,mRef,mass,massResol);
2238  // P = GL->Integral(mass-5*massResol, mass+5*massResol);
2239 
2240  Int_t np = 100;
2241  double *x = new double[np];
2242  double *w = new double[np];
2243  GL->SetParameters(ResGamma[ires], ResMass[ires], mass, massResol);
2244  GL->CalcGaussLegendreSamplingPoints(np, x, w, 0.1e-15);
2245  P = GL->IntegralFast(np, x, w, ResMass[ires] - 10 * ResGamma[ires], ResMass[ires] + 10 * ResGamma[ires]);
2246  delete[] x;
2247  delete[] w;
2248 
2249  // If we are too far away we set P to epsilon and forget about this event
2250  // ----------------------------------------------------------------------
2251  if (P < 1.0e-12) {
2252  P = 1.0e-12;
2253  if (debug > 9)
2254  std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2255  << massResol << ": used epsilon" << std::endl;
2256  return P;
2257  }
2258 
2259  if (debug > 9)
2260  std::cout << "Mass, gamma , mref, width, P: " << mass << " " << ResGamma[ires] << " " << ResMass[ires] << " "
2261  << massResol << " " << P << std::endl;
2262  return P;
2263 }
const double Pi
T w() const
static double x[7][10000]
static int debug
static double ResMass[6]
TF1 * GL
int np
Definition: AMPTWrapper.h:43
static double ResGamma[6]
std::pair< OmniClusterRef, TrackingParticleRef > P

◆ massProb() [2/3]

double MuScleFitUtils::massProb ( const double &  mass,
const double &  resEta,
const double &  rapidity,
const double &  massResol,
const std::vector< double > &  parval,
const bool  doUseBkgrWindow,
const double &  eta1,
const double &  eta2 
)
static

Definition at line 710 of file MuScleFitUtils.cc.

References HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, EcalPhiSymFlatTableProducers_cfi::id, createfilelist::int, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, EgHLTOffHistBins_cfi::mass, massProb(), and AlCaHLTBitMon_ParallelJobs::p.

717  {
718 #ifdef USE_CALLGRIND
719  CALLGRIND_START_INSTRUMENTATION;
720 #endif
721 
722  double *p = new double[(int)(parval.size())];
723  // Replaced by auto_ptr, which handles delete at the end
724  // Removed auto_ptr, check massResolution for an explanation.
725  // std::unique_ptr<double> p(new double[(int)(parval.size())]);
726 
727  std::vector<double>::const_iterator it = parval.begin();
728  int id = 0;
729  for (; it != parval.end(); ++it, ++id) {
730  // (&*p)[id] = *it;
731  p[id] = *it;
732  }
733  // p must be passed by value as below:
734  double massProbability = massProb(mass, resEta, rapidity, massResol, p, doUseBkgrWindow, eta1, eta2);
735  delete[] p;
736 
737 #ifdef USE_CALLGRIND
738  CALLGRIND_STOP_INSTRUMENTATION;
739  CALLGRIND_DUMP_STATS;
740 #endif
741 
742  return massProbability;
743 }
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)

◆ massProb() [3/3]

double MuScleFitUtils::massProb ( const double &  mass,
const double &  resEta,
const double &  rapidity,
const double &  massResol,
double *  parval,
const bool  doUseBkgrWindow,
const double &  eta1,
const double &  eta2 
)
static

Definition at line 889 of file MuScleFitUtils.cc.

References funct::abs(), BackgroundHandler::backgroundFunction(), backgroundHandler, backgroundProb_, checkMassWindow(), gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, GLNorm, GLValue, GLZNorm, GLZValue, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, loopCounter, EgHLTOffHistBins_cfi::mass, minuitLoop_, MuonType, parBgr, CrossSectionHandler::parNum(), parResol, parScale, probability(), rapidityBinsForZ_, CrossSectionHandler::relativeCrossSections(), resfind, ResHalfWidth, BackgroundHandler::resMass(), ResMass, signalProb_, totalResNum, and BackgroundHandler::windowBorders().

896  {
897  // This routine computes the likelihood that a given measured mass "measMass" is
898  // the result of a reference mass ResMass[] if the resolution
899  // expected for the two muons is massResol.
900  // This version includes two parameters (the last two in parval, by default)
901  // to size up the background fraction and its relative normalization with respect
902  // to the signal shape.
903  //
904  // We model the signal probability with a Lorentz L(M,H) of resonance mass M and natural width H
905  // convoluted with a gaussian G(m,s) of measured mass m and expected mass resolution s,
906  // by integrating over the intersection of the supports of L and G (which can be made to coincide with
907  // the region where L is non-zero, given that H<<s in most cases) the product L(M,H)*G(m,s) dx as follows:
908  //
909  // GL(m,s) = Int(M-10H,M+10H) [ L(x-M,H) * G(x-m,s) ] dx
910  //
911  // The above convolution is computed numerically by an independent root macro, Probs.C, which outputs
912  // the values in six 1001x1001 grids, one per resonance.
913  //
914  // NB THe following block of explanations for background models is outdated, see detailed
915  // explanations where the code computes PB.
916  // +++++++++++++++++++++++
917  // For the background, instead, we have two choices: a linear and an exponential model.
918  // * For the linear model, we choose a one-parameter form whereby the line is automatically normalized
919  // in the support [x1,x2] where we defined our "signal region", as follows:
920  //
921  // B(x;b) = 1/(x2-x1) + {x - (x2+x1)/2} * b
922  //
923  // Defined as above, B(x) is a line passing for the point of coordinates (x1+x2)/2, 1/(x2-x1),
924  // whose slope b has as support the interval ( -2/[(x1-x2)*(x1+x2)], 2/[(x1-x2)*(x1+x2)] )
925  // so that B(x) is always positive-definite in [x1,x2].
926  //
927  // * For the exponential model, we define B(x;b) as
928  //
929  // B(x;b) = b * { exp(-b*x) / [exp(-b*x1)-exp(-b*x2)] }
930  //
931  // This one-parameter definition is automatically normalized to unity in [x1,x2], with a parameter
932  // b which has to be positive in order for the slope to be negative.
933  // Please note that this model is not useful in most circumstances; a more useful form would be one
934  // which included a linear component.
935  // ++++++++++++++++++++++
936  //
937  // Once GL(m,s) and B(x;b) are defined, we introduce a further parameter a, such that we can have the
938  // likelihood control the relative fraction of signal and background. We first normalize GL(m,s) for
939  // any given s by taking the integral
940  //
941  // Int(x1,x2) GL(m,s) dm = K_s
942  //
943  // We then define the probability as
944  //
945  // P(m,s,a,b) = GL(m,s)/K_s * a + B(x,b) * (1-a)
946  //
947  // with a taking on values in the interval [0,1].
948  // Defined as above, the probability is well-behaved, in the sense that it has a value between 0 and 1,
949  // and the four parameters m,s,a,b fully control its shape.
950  //
951  // It is to be noted that the formulation above requires the computation of two rather time-consuming
952  // integrals. The one defining GL(m,s) can be stored in a TH2D and loaded by the constructor from a
953  // file suitably prepared, and this will save loads of computing time.
954  // ----------------------------------------------------------------------------------------------------
955 
956  double P = 0.;
957  int crossSectionParShift = parResol.size() + parScale.size();
958  // Take the relative cross sections
959  std::vector<double> relativeCrossSections =
960  crossSectionHandler->relativeCrossSections(&(parval[crossSectionParShift]), resfind);
961  // for( unsigned int i=0; i<relativeCrossSections.size(); ++i ) {
962  // std::cout << "relativeCrossSections["<<i<<"] = " << relativeCrossSections[i] << std::endl;
963  // std::cout << "parval["<<crossSectionParShift+i<<"] = " << parval[crossSectionParShift+i] << std::endl;
964  // }
965 
966  // int bgrParShift = crossSectionParShift + parCrossSection.size();
967  int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
968  double Bgrp1 = 0.;
969  // double Bgrp2 = 0.;
970  // double Bgrp3 = 0.;
971 
972  // NB defined as below, P is a non-rigorous "probability" that we observe
973  // a dimuon mass "mass", given ResMass[], gamma, and massResol. It is what we need for the
974  // fit which finds the best resolution parameters, though. A definition which is
975  // more properly a probability is given below (in massProb2()), where the result
976  // cannot be used to fit resolution parameters because the fit would always prefer
977  // to set the res parameters to the minimum possible value (best resolution),
978  // to have a probability close to one of observing any mass.
979  // -------------------------------------------------------------------------------
980 
981  // Determine what resonance(s) we have to deal with
982  // NB for now we assume equal xs for each resonance
983  // so we do not assign them different weights
984  // ------------------------------------------------
985  double PS[6] = {0.};
986  double PB = 0.;
987  double PStot[6] = {0.};
988 
989  // Should be removed because it is not used
990  bool resConsidered[6] = {false};
991 
992  bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
993  // bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
994 
995  // First check the Z, which is divided in 24 rapidity bins
996  // NB max value of Z rapidity to be considered is 2.4 here
997  // -------------------------------------------------------
998 
999  // Do this only if we want to use the rapidity bins for the Z
1001  // ATTENTION: cut on Z rapidity at 2.4 since we only have histograms up to that value
1002  // std::pair<double, double> windowFactors = backgroundHandler->windowFactors( useBackgroundWindow, 0 );
1003  std::pair<double, double> windowBorders = backgroundHandler->windowBorders(useBackgroundWindow, 0);
1004  if (resfind[0] > 0
1005  // && checkMassWindow( mass, 0,
1006  // backgroundHandler->resMass( useBackgroundWindow, 0 ),
1007  // windowFactors.first, windowFactors.second )
1008  && checkMassWindow(mass, windowBorders.first, windowBorders.second)
1009  // && std::abs(rapidity)<2.4
1010  ) {
1011  int iY = (int)(std::abs(rapidity) * 10.);
1012  if (iY > 23)
1013  iY = 23;
1014 
1015  if (MuScleFitUtils::debug > 1)
1016  std::cout << "massProb:resFound = 0, rapidity bin =" << iY << std::endl;
1017 
1018  // In this case the last value is the rapidity bin
1019  PS[0] = probability(mass, massResol, GLZValue, GLZNorm, 0, iY);
1020 
1021  if (PS[0] != PS[0]) {
1022  std::cout << "ERROR: PS[0] = nan, setting it to 0" << std::endl;
1023  PS[0] = 0;
1024  }
1025 
1026  // std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
1027  // &(parval[bgrParShift]), MuScleFitUtils::totalResNum, 0,
1028  // resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
1029 
1030  std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction(doBackgroundFit[loopCounter],
1031  &(parval[bgrParShift]),
1033  0,
1034  resConsidered,
1035  ResMass,
1036  ResHalfWidth,
1037  MuonType,
1038  mass,
1039  eta1,
1040  eta2);
1041 
1042  Bgrp1 = bgrResult.first;
1043  // When fitting the background we have only one Bgrp1
1044  // When not fitting the background we have many only in a superposition region and this case is treated
1045  // separately after this loop
1046  PB = bgrResult.second;
1047  if (PB != PB)
1048  PB = 0;
1049  PStot[0] = (1 - Bgrp1) * PS[0] + Bgrp1 * PB;
1050 
1051  // PStot[0] *= crossSectionHandler->crossSection(0);
1052  // PStot[0] *= parval[crossSectionParShift];
1053  PStot[0] *= relativeCrossSections[0];
1054  if (MuScleFitUtils::debug > 0)
1055  std::cout << "PStot[" << 0 << "] = "
1056  << "(1-" << Bgrp1 << ")*" << PS[0] << " + " << Bgrp1 << "*" << PB << " = " << PStot[0]
1057  << " (reletiveXS = )" << relativeCrossSections[0] << std::endl;
1058  } else {
1059  if (debug > 0) {
1060  std::cout << "Mass = " << mass << " outside range with rapidity = " << rapidity << std::endl;
1061  std::cout << "with resMass = " << backgroundHandler->resMass(useBackgroundWindow, 0)
1062  << " and left border = " << windowBorders.first << " right border = " << windowBorders.second
1063  << std::endl;
1064  }
1065  }
1066  }
1067  // Next check the other resonances
1068  // -------------------------------
1069  int firstRes = 1;
1071  firstRes = 0;
1072  for (int ires = firstRes; ires < 6; ++ires) {
1073  if (resfind[ires] > 0) {
1074  // First is left, second is right (returns (1,1) in the case of resonances, it could be improved avoiding the call in this case)
1075  // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
1076  std::pair<double, double> windowBorder = backgroundHandler->windowBorders(useBackgroundWindow, ires);
1077  if (checkMassWindow(mass, windowBorder.first, windowBorder.second)) {
1078  if (MuScleFitUtils::debug > 1)
1079  std::cout << "massProb:resFound = " << ires << std::endl;
1080 
1081  // In this case the rapidity value is instead the resonance index again.
1082  PS[ires] = probability(mass, massResol, GLValue, GLNorm, ires, ires);
1083 
1084  std::pair<double, double> bgrResult =
1086  &(parval[bgrParShift]),
1088  ires,
1089  // resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
1090  resConsidered,
1091  ResMass,
1092  ResHalfWidth,
1093  MuonType,
1094  mass,
1095  eta1,
1096  eta2);
1097  Bgrp1 = bgrResult.first;
1098  PB = bgrResult.second;
1099 
1100  if (PB != PB)
1101  PB = 0;
1102  PStot[ires] = (1 - Bgrp1) * PS[ires] + Bgrp1 * PB;
1103  if (MuScleFitUtils::debug > 0)
1104  std::cout << "PStot[" << ires << "] = "
1105  << "(1-" << Bgrp1 << ")*" << PS[ires] << " + " << Bgrp1 << "*" << PB << " = " << PStot[ires]
1106  << std::endl;
1107 
1108  PStot[ires] *= relativeCrossSections[ires];
1109  }
1110  }
1111  }
1112 
1113  for (int i = 0; i < 6; ++i) {
1114  P += PStot[i];
1115  }
1116 
1117  if (MuScleFitUtils::signalProb_ != nullptr && MuScleFitUtils::backgroundProb_ != nullptr) {
1118  double PStotTemp = 0.;
1119  for (int i = 0; i < 6; ++i) {
1120  PStotTemp += PS[i] * relativeCrossSections[i];
1121  }
1122  if (PStotTemp != PStotTemp) {
1123  std::cout << "ERROR: PStotTemp = nan!!!!!!!!!" << std::endl;
1124  int parnumber = (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size());
1125  for (int i = 0; i < 6; ++i) {
1126  std::cout << "PS[i] = " << PS[i] << std::endl;
1127  if (PS[i] != PS[i]) {
1128  std::cout << "mass = " << mass << std::endl;
1129  std::cout << "massResol = " << massResol << std::endl;
1130  for (int j = 0; j < parnumber; ++j) {
1131  std::cout << "parval[" << j << "] = " << parval[j] << std::endl;
1132  }
1133  }
1134  }
1135  }
1136  if (PStotTemp == PStotTemp) {
1137  MuScleFitUtils::signalProb_->SetBinContent(
1139  MuScleFitUtils::signalProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PStotTemp);
1140  }
1141  if (debug > 0)
1142  std::cout << "mass = " << mass << ", P = " << P << ", PStot = " << PStotTemp << ", PB = " << PB
1143  << ", bgrp1 = " << Bgrp1 << std::endl;
1144 
1145  MuScleFitUtils::backgroundProb_->SetBinContent(
1147  }
1148  return P;
1149 }
static double GLValue[6][1001][1001]
static unsigned int loopCounter
static std::vector< int > doBackgroundFit
static std::vector< double > parResol
static TH1D * backgroundProb_
static BackgroundHandler * backgroundHandler
std::pair< double, double > backgroundFunction(const bool doBackgroundFit, const double *parval, const int resTotNum, const int ires, const bool *resConsidered, const double *ResMass, const double ResHalfWidth[], const int MuonType, const double &mass, const double &eta1, const double &eta2)
static int debug
static double ResMass[6]
static bool checkMassWindow(const double &mass, const double &leftBorder, const double &rightBorder)
Method to check if the mass value is within the mass window of the i-th resonance.
static const int totalResNum
static int minuitLoop_
static double GLZNorm[40][1001]
static std::vector< double > parBgr
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
static TH1D * signalProb_
static double GLZValue[40][1001][1001]
double resMass(const bool doBackgroundFit, const int ires)
static std::vector< double > parScale
std::vector< double > relativeCrossSections(const double *variables, const std::vector< int > &resfind)
Perform a variable transformation from N-1 to relative cross sections.
std::pair< OmniClusterRef, TrackingParticleRef > P
static bool rapidityBinsForZ_
static double GLNorm[6][1001]
std::pair< double, double > windowBorders(const bool doBackgroundFit, const int ires)
Returns the appropriate window borders depending on whether the background is being fitted and on the...
static int MuonType
static double probability(const double &mass, const double &massResol, const double GLvalue[][1001][1001], const double GLnorm[][1001], const int iRes, const int iY)
Computes the probability given the mass, mass resolution and the arrays with the probabilities and th...
static std::vector< int > resfind
static double ResHalfWidth[6]
static CrossSectionHandler * crossSectionHandler

◆ massResolution() [1/5]

static double MuScleFitUtils::massResolution ( const lorentzVector mu1,
const lorentzVector mu2 
)
static

◆ massResolution() [2/5]

double MuScleFitUtils::massResolution ( const lorentzVector mu1,
const lorentzVector mu2,
const std::vector< double > &  parval 
)
static

Definition at line 524 of file MuScleFitUtils.cc.

References EcalPhiSymFlatTableProducers_cfi::id, createfilelist::int, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, massResolution(), and AlCaHLTBitMon_ParallelJobs::p.

526  {
527  // double * p = new double[(int)(parval.size())];
528  // Replaced by auto_ptr, which handles delete at the end
529  // --------- //
530  // ATTENTION //
531  // --------- //
532  // auto_ptr calls delete, not delete[] and thus it must
533  // not be used with arrays. There are alternatives see
534  // e.g.: http://www.gotw.ca/gotw/042.htm. The best
535  // alternative seems to be to switch to vector though.
536  // std::unique_ptr<double> p(new double[(int)(parval.size())]);
537 
538  double *p = new double[(int)(parval.size())];
539  std::vector<double>::const_iterator it = parval.begin();
540  int id = 0;
541  for (; it != parval.end(); ++it, ++id) {
542  // (&*p)[id] = *it;
543  p[id] = *it;
544  }
545  double massRes = massResolution(mu1, mu2, p);
546  delete[] p;
547  return massRes;
548 }
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)

◆ massResolution() [3/5]

static double MuScleFitUtils::massResolution ( const lorentzVector mu1,
const lorentzVector mu2,
std::unique_ptr< double >  parval 
)
static

◆ massResolution() [4/5]

double MuScleFitUtils::massResolution ( const lorentzVector mu1,
const lorentzVector mu2,
double *  parval 
)
static

We use the following formula:

M = sqrt ( (E1+E2)^2 - (P1+P2)^2 )

where we express E and P as a function of Pt, phi, and theta:

E = sqrt ( Pt^2*(1+cotg(theta)^2) + M_mu^2 )
Px = Pt*cos(phi), Py = Pt*sin(phi), Pz = Pt*cotg(theta)

from which we find

M = sqrt( 2*M_mu^2 + 2*sqrt(Pt1^2/sin(theta1)^2 + M_mu^2)*sqrt(Pt2^2/sin(theta2)^2 + M_mu^2) - 2*Pt1*Pt2* ( cos(phi1-phi2) + cotg(theta1)*cotg(theta2) ) )

and derive WRT Pt1, Pt2, phi1, phi2, theta1, theta2 to get the resolution.

Definition at line 567 of file MuScleFitUtils.cc.

References funct::abs(), funct::cos(), counter_resprob, gather_cfg::cout, resolutionFunctionBase< T >::covPt1Pt2(), debug, debugMassResol_, MuScleFitUtils::massResolComponentsStruct::dmdcotgth1, MuScleFitUtils::massResolComponentsStruct::dmdcotgth2, MuScleFitUtils::massResolComponentsStruct::dmdphi1, MuScleFitUtils::massResolComponentsStruct::dmdphi2, MuScleFitUtils::massResolComponentsStruct::dmdpt1, MuScleFitUtils::massResolComponentsStruct::dmdpt2, HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, JetChargeProducer_cfi::exp, LogDebug, EgHLTOffHistBins_cfi::mass, massResolComponents, mMu2, funct::pow(), HLT_2024v14_cff::pt1, HLT_2024v14_cff::pt2, resfind, ResHalfWidth, ResMass, ResMaxSigma, resolutionFunction, resolutionFunctionBase< T >::sigmaCotgTh(), resolutionFunctionBase< T >::sigmaPhi(), resolutionFunctionBase< T >::sigmaPt(), funct::sin(), and mathSSE::sqrt().

567  {
568  double mass = (mu1 + mu2).mass();
569  double pt1 = mu1.Pt();
570  double phi1 = mu1.Phi();
571  double eta1 = mu1.Eta();
572  double theta1 = 2 * atan(exp(-eta1));
573  double pt2 = mu2.Pt();
574  double phi2 = mu2.Phi();
575  double eta2 = mu2.Eta();
576  double theta2 = 2 * atan(exp(-eta2));
577 
578  double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
579  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
580  pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
581  mass;
582  double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
583  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
584  pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
585  mass;
586  double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
587  double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
588  double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
589  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
590  pt1 * pt2 * cos(theta2) / sin(theta2)) /
591  mass;
592  double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
593  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
594  pt2 * pt1 * cos(theta1) / sin(theta1)) /
595  mass;
596 
597  if (debugMassResol_) {
598  massResolComponents.dmdpt1 = dmdpt1;
599  massResolComponents.dmdpt2 = dmdpt2;
600  massResolComponents.dmdphi1 = dmdphi1;
601  massResolComponents.dmdphi2 = dmdphi2;
602  massResolComponents.dmdcotgth1 = dmdcotgth1;
603  massResolComponents.dmdcotgth2 = dmdcotgth2;
604  }
605 
606  // Resolution parameters:
607  // ----------------------
608  double sigma_pt1 = resolutionFunction->sigmaPt(pt1, eta1, parval);
609  double sigma_pt2 = resolutionFunction->sigmaPt(pt2, eta2, parval);
610  double sigma_phi1 = resolutionFunction->sigmaPhi(pt1, eta1, parval);
611  double sigma_phi2 = resolutionFunction->sigmaPhi(pt2, eta2, parval);
612  double sigma_cotgth1 = resolutionFunction->sigmaCotgTh(pt1, eta1, parval);
613  double sigma_cotgth2 = resolutionFunction->sigmaCotgTh(pt2, eta2, parval);
614  double cov_pt1pt2 = resolutionFunction->covPt1Pt2(pt1, eta1, pt2, eta2, parval);
615 
616  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
617  // multiply it by pt.
618  double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
619  std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
620  std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2) +
621  2 * dmdpt1 * dmdpt2 * cov_pt1pt2 * sigma_pt1 * sigma_pt2);
622 
623  if (debug > 19) {
624  std::cout << " Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1) / sin(theta1) << " - Pt2=" << pt2
625  << " phi2=" << phi2 << " cotgth2=" << cos(theta2) / sin(theta2) << std::endl;
626  std::cout << " P[0]=" << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3]
627  << std::endl;
628  std::cout << " Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1=" << sigma_pt1
629  << " sigma_pt2=" << sigma_pt2 << std::endl;
630  std::cout << " Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1=" << sigma_phi1
631  << " sigma_phi2=" << sigma_phi2 << std::endl;
632  std::cout << " Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2 << " sigma_cotgth1=" << sigma_cotgth1
633  << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
634  std::cout << " Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2 << " : " << mass << " +- "
635  << mass_res << std::endl;
636  }
637 
638  // Debug std::cout
639  // ----------
640  bool didit = false;
641  for (int ires = 0; ires < 6; ires++) {
642  if (!didit && resfind[ires] > 0 && std::abs(mass - ResMass[ires]) < ResHalfWidth[ires]) {
643  if (mass_res > ResMaxSigma[ires] && counter_resprob < 100) {
644  counter_resprob++;
645  LogDebug("MuScleFitUtils") << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
646  didit = true;
647  }
648  }
649  }
650 
651  return mass_res;
652 }
static bool debugMassResol_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
static int debug
static double ResMass[6]
static struct MuScleFitUtils::massResolComponentsStruct massResolComponents
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual double covPt1Pt2(const double &pt1, const double &eta1, const double &pt2, const double &eta2, const T &parval)
Definition: Functions.h:834
static double ResMaxSigma[6]
virtual double sigmaCotgTh(const double &pt, const double &eta, const T &parval)=0
static const double mMu2
virtual double sigmaPhi(const double &pt, const double &eta, const T &parval)=0
static std::vector< int > resfind
static double ResHalfWidth[6]
static int counter_resprob
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
static resolutionFunctionBase< double * > * resolutionFunction
#define LogDebug(id)

◆ massResolution() [5/5]

double MuScleFitUtils::massResolution ( const lorentzVector mu1,
const lorentzVector mu2,
const ResolutionFunction resolFunc 
)
static

This method can be used outside MuScleFit. It gets the ResolutionFunction that must have been built with the parameters.
TO-DO: this method duplicates the code in the previous method. It should be changed to avoid the duplication.

Definition at line 658 of file MuScleFitUtils.cc.

References funct::cos(), HLT_2024v14_cff::eta1, HLT_2024v14_cff::eta2, JetChargeProducer_cfi::exp, EgHLTOffHistBins_cfi::mass, mMu2, funct::pow(), HLT_2024v14_cff::pt1, HLT_2024v14_cff::pt2, ResolutionFunction::sigmaCotgTh(), ResolutionFunction::sigmaPhi(), ResolutionFunction::sigmaPt(), funct::sin(), and mathSSE::sqrt().

660  {
661  double mass = (mu1 + mu2).mass();
662  double pt1 = mu1.Pt();
663  double phi1 = mu1.Phi();
664  double eta1 = mu1.Eta();
665  double theta1 = 2 * atan(exp(-eta1));
666  double pt2 = mu2.Pt();
667  double phi2 = mu2.Phi();
668  double eta2 = mu2.Eta();
669  double theta2 = 2 * atan(exp(-eta2));
670 
671  double dmdpt1 = (pt1 / std::pow(sin(theta1), 2) *
672  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
673  pt2 * (cos(phi1 - phi2) + cos(theta1) * cos(theta2) / (sin(theta1) * sin(theta2)))) /
674  mass;
675  double dmdpt2 = (pt2 / std::pow(sin(theta2), 2) *
676  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
677  pt1 * (cos(phi2 - phi1) + cos(theta2) * cos(theta1) / (sin(theta2) * sin(theta1)))) /
678  mass;
679  double dmdphi1 = pt1 * pt2 / mass * sin(phi1 - phi2);
680  double dmdphi2 = pt2 * pt1 / mass * sin(phi2 - phi1);
681  double dmdcotgth1 = (pt1 * pt1 * cos(theta1) / sin(theta1) *
682  sqrt((std::pow(pt2 / sin(theta2), 2) + mMu2) / (std::pow(pt1 / sin(theta1), 2) + mMu2)) -
683  pt1 * pt2 * cos(theta2) / sin(theta2)) /
684  mass;
685  double dmdcotgth2 = (pt2 * pt2 * cos(theta2) / sin(theta2) *
686  sqrt((std::pow(pt1 / sin(theta1), 2) + mMu2) / (std::pow(pt2 / sin(theta2), 2) + mMu2)) -
687  pt2 * pt1 * cos(theta1) / sin(theta1)) /
688  mass;
689 
690  // Resolution parameters:
691  // ----------------------
692  double sigma_pt1 = resolFunc.sigmaPt(mu1);
693  double sigma_pt2 = resolFunc.sigmaPt(mu2);
694  double sigma_phi1 = resolFunc.sigmaPhi(mu1);
695  double sigma_phi2 = resolFunc.sigmaPhi(mu2);
696  double sigma_cotgth1 = resolFunc.sigmaCotgTh(mu1);
697  double sigma_cotgth2 = resolFunc.sigmaCotgTh(mu2);
698 
699  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
700  // multiply it by pt.
701  double mass_res = sqrt(std::pow(dmdpt1 * sigma_pt1 * pt1, 2) + std::pow(dmdpt2 * sigma_pt2 * pt2, 2) +
702  std::pow(dmdphi1 * sigma_phi1, 2) + std::pow(dmdphi2 * sigma_phi2, 2) +
703  std::pow(dmdcotgth1 * sigma_cotgth1, 2) + std::pow(dmdcotgth2 * sigma_cotgth2, 2));
704 
705  return mass_res;
706 }
double sigmaCotgTh(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double sigmaPt(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
static const double mMu2
double sigmaPhi(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ minimizeLikelihood()

void MuScleFitUtils::minimizeLikelihood ( )
static

Definition at line 1197 of file MuScleFitUtils.cc.

References funct::abs(), backgroundHandler, backgroundProb_, BgrFitType, svgfig::canvas(), RPCNoise_example::check, checkMassWindow(), computeMinosErrors_, gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, doCrossSectionFit, doResolFit, doScaleFit, duringMinos_, FitStrategy, mps_fire::i, createfilelist::int, invDimuonMass(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, likelihood(), likelihoodInLoop_, loopCounter, EgHLTOffHistBins_cfi::mass, massWindowHalfWidth, minimumShapePlots_, minuitLoop_, MuonType, mergeVDriftHistosByStation::name, makeMEIFBenchmarkPlots::nev, normalizationChanged_, parBgr, parBgrFix, parBgrOrder, parCrossSection, parCrossSectionOrder, parfix, scaleFunctionBase< T >::parNum(), CrossSectionHandler::parNum(), resolutionFunctionBase< T >::parNum(), parorder, parResol, parResolFix, parResolMax, parResolMin, parResolOrder, parResolStep, parScale, parScaleFix, parScaleMax, parScaleMin, parScaleOrder, parScaleStep, parvalue, TrackingDataMCValidation_Standalone_cff::pmin, ReducedSavedPair, BackgroundHandler::regionsParNum(), CrossSectionHandler::relativeCrossSections(), CrossSectionHandler::releaseParameters(), BackgroundHandler::rescale(), scaleFunctionBase< T >::resetParameters(), resfind, ResMass, ResolFitType, resolutionFunctionForVec, rminPtr_, SavedPair, scaleFitNotDone_, ScaleFitType, scaleFunction, scaleFunctionForVec, scaleFunctionBase< T >::setParameters(), BackgroundHandler::setParameters(), CrossSectionHandler::setParameters(), resolutionFunctionBase< T >::setParameters(), signalProb_, startWithSimplex_, BackgroundHandler::unlockParameter(), and BackgroundHandler::windowBorders().

Referenced by MuScleFit::endOfFastLoop().

1197  {
1198  // Output file with fit parameters resulting from minimization
1199  // -----------------------------------------------------------
1200  std::ofstream FitParametersFile;
1201  FitParametersFile.open("FitParameters.txt", std::ios::app);
1202  FitParametersFile << "Fitting with resolution, scale, bgr function # " << ResolFitType << " " << ScaleFitType << " "
1203  << BgrFitType << " - Iteration " << loopCounter << std::endl;
1204 
1205  // Fill parvalue and other vectors needed for the fitting
1206  // ------------------------------------------------------
1207 
1208  // ----- //
1209  // FIXME //
1210  // ----- //
1211  // this was changed to verify the possibility that fixed parameters influence the errors.
1212  // It must be 0 otherwise the parameters for resonances will not be passed by minuit (will be always 0).
1213  // Should be removed.
1214  int parForResonanceWindows = 0;
1215  // int parnumber = (int)(parResol.size()+parScale.size()+parCrossSection.size()+parBgr.size() - parForResonanceWindows);
1216  int parnumber =
1217  (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size() - parForResonanceWindows);
1218 
1219  int parnumberAll = (int)(parResol.size() + parScale.size() + crossSectionHandler->parNum() + parBgr.size());
1220 
1221  // parvalue is a std::vector<std::vector<double> > storing all the parameters from all the loops
1222  parvalue.push_back(parResol);
1223  std::vector<double> *tmpVec = &(parvalue.back());
1224 
1225  // If this is not the first loop we want to start from neutral values
1226  // Otherwise the scale will start with values correcting again a bias
1227  // that is already corrected.
1228  if (scaleFitNotDone_) {
1229  tmpVec->insert(tmpVec->end(), parScale.begin(), parScale.end());
1230  std::cout << "scaleFitNotDone: tmpVec->size() = " << tmpVec->size() << std::endl;
1231  } else {
1232  scaleFunction->resetParameters(tmpVec);
1233  std::cout << "scaleFitDone: tmpVec->size() = " << tmpVec->size() << std::endl;
1234  }
1235  tmpVec->insert(tmpVec->end(), parCrossSection.begin(), parCrossSection.end());
1236  tmpVec->insert(tmpVec->end(), parBgr.begin(), parBgr.end());
1237  int i = 0;
1238  std::vector<double>::const_iterator it = tmpVec->begin();
1239  for (; it != tmpVec->end(); ++it, ++i) {
1240  std::cout << "tmpVec[" << i << "] = " << *it << std::endl;
1241  }
1242 
1243  // Empty vector of size = number of cross section fitted parameters. Note that the cross section
1244  // fit works in a different way than the others and it uses ratios of the paramters passed via cfg.
1245  // We use this empty vector for compatibility with the rest of the structure.
1246  std::vector<int> crossSectionParNumSizeVec(MuScleFitUtils::crossSectionHandler->parNum(), 0);
1247 
1248  std::vector<int> parfix(parResolFix);
1249  parfix.insert(parfix.end(), parScaleFix.begin(), parScaleFix.end());
1250  parfix.insert(parfix.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end());
1251  parfix.insert(parfix.end(), parBgrFix.begin(), parBgrFix.end());
1252 
1253  std::vector<int> parorder(parResolOrder);
1254  parorder.insert(parorder.end(), parScaleOrder.begin(), parScaleOrder.end());
1255  parorder.insert(parorder.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end());
1256  parorder.insert(parorder.end(), parBgrOrder.begin(), parBgrOrder.end());
1257 
1258  // This is filled later
1259  std::vector<double> parerr(3 * parnumberAll, 0.);
1260 
1261  if (debug > 19) {
1262  std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters before likelihood " << std::endl;
1263  for (unsigned int i = 0; i < (unsigned int)parnumberAll; i++) {
1264  std::cout << " Par # " << i << " = " << parvalue[loopCounter][i] << " : free = " << parfix[i]
1265  << "; order = " << parorder[i] << std::endl;
1266  }
1267  }
1268 
1269  // // Background rescaling from regions to resonances
1270  // // -----------------------------------------------
1271  // // If we are in a loop > 0 and we are not fitting the background, but we have fitted it in the previous iteration
1272  // if( loopCounter > 0 && !(doBackgroundFit[loopCounter]) && doBackgroundFit[loopCounter-1] ) {
1273  // // This rescales from regions to resonances
1274  // int localMuonType = MuonType;
1275  // if( MuonType > 2 ) localMuonType = 2;
1276  // backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
1277  // MuScleFitUtils::SavedPair);
1278  // }
1279 
1280  // Init Minuit
1281  // -----------
1282  TMinuit rmin(parnumber);
1283  rminPtr_ = &rmin;
1284  rmin.SetFCN(likelihood); // Unbinned likelihood
1285  // Standard initialization of minuit parameters:
1286  // sets input to be $stdin, output to be $stdout
1287  // and saving to a file.
1288  rmin.mninit(5, 6, 7);
1289  int ierror = 0;
1290  int istat;
1291  double arglis[4];
1292  arglis[0] = FitStrategy; // Strategy 1 or 2
1293  // 1 standard
1294  // 2 try to improve minimum (slower)
1295  rmin.mnexcm("SET STR", arglis, 1, ierror);
1296 
1297  arglis[0] = 10001;
1298  // Set the random seed for the generator used in SEEk to a fixed value for reproducibility
1299  rmin.mnexcm("SET RAN", arglis, 1, ierror);
1300 
1301  // Set fit parameters
1302  // ------------------
1303  double *Start = new double[parnumberAll];
1304  double *Step = new double[parnumberAll];
1305  double *Mini = new double[parnumberAll];
1306  double *Maxi = new double[parnumberAll];
1307  int *ind = new int[parnumberAll]; // Order of release of parameters
1308  TString *parname = new TString[parnumberAll];
1309 
1310  if (!parResolStep.empty() && !parResolMin.empty() && !parResolMax.empty()) {
1312  Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parResolStep, parResolMin, parResolMax, MuonType);
1313  } else {
1315  Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, MuonType);
1316  }
1317 
1318  // Take the number of parameters in the resolutionFunction and displace the arrays passed to the scaleFunction
1320 
1321  if (!parScaleStep.empty() && !parScaleMin.empty() && !parScaleMax.empty()) {
1323  &(Step[resParNum]),
1324  &(Mini[resParNum]),
1325  &(Maxi[resParNum]),
1326  &(ind[resParNum]),
1327  &(parname[resParNum]),
1328  parScale,
1329  parScaleOrder,
1330  parScaleStep,
1331  parScaleMin,
1332  parScaleMax,
1333  MuonType);
1334  } else {
1336  &(Step[resParNum]),
1337  &(Mini[resParNum]),
1338  &(Maxi[resParNum]),
1339  &(ind[resParNum]),
1340  &(parname[resParNum]),
1341  parScale,
1342  parScaleOrder,
1343  MuonType);
1344  }
1345 
1346  // Initialize cross section parameters
1347  int crossSectionParShift = resParNum + MuScleFitUtils::scaleFunctionForVec->parNum();
1348  MuScleFitUtils::crossSectionHandler->setParameters(&(Start[crossSectionParShift]),
1349  &(Step[crossSectionParShift]),
1350  &(Mini[crossSectionParShift]),
1351  &(Maxi[crossSectionParShift]),
1352  &(ind[crossSectionParShift]),
1353  &(parname[crossSectionParShift]),
1356  resfind);
1357 
1358  // Initialize background parameters
1359  int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
1360  MuScleFitUtils::backgroundHandler->setParameters(&(Start[bgrParShift]),
1361  &(Step[bgrParShift]),
1362  &(Mini[bgrParShift]),
1363  &(Maxi[bgrParShift]),
1364  &(ind[bgrParShift]),
1365  &(parname[bgrParShift]),
1366  parBgr,
1367  parBgrOrder,
1368  MuonType);
1369 
1370  for (int ipar = 0; ipar < parnumber; ++ipar) {
1371  std::cout << "parname[" << ipar << "] = " << parname[ipar] << std::endl;
1372  std::cout << "Start[" << ipar << "] = " << Start[ipar] << std::endl;
1373  std::cout << "Step[" << ipar << "] = " << Step[ipar] << std::endl;
1374  std::cout << "Mini[" << ipar << "] = " << Mini[ipar] << std::endl;
1375  std::cout << "Maxi[" << ipar << "] = " << Maxi[ipar] << std::endl;
1376 
1377  rmin.mnparm(ipar, parname[ipar], Start[ipar], Step[ipar], Mini[ipar], Maxi[ipar], ierror);
1378 
1379  // Testing without limits
1380  // rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], 0, 0, ierror );
1381  }
1382 
1383  // Do minimization
1384  // ---------------
1385  if (debug > 19)
1386  std::cout << "[MuScleFitUtils-minimizeLikelihood]: Starting minimization" << std::endl;
1387  double fmin;
1388  double fdem;
1389  double errdef;
1390  int npari;
1391  int nparx;
1392  rmin.mnexcm("CALL FCN", arglis, 1, ierror);
1393 
1394  // First, fix all parameters
1395  // -------------------------
1396  if (debug > 19)
1397  std::cout << "[MuScleFitUtils-minimizeLikelihood]: First fix all parameters ...";
1398  for (int ipar = 0; ipar < parnumber; ipar++) {
1399  rmin.FixParameter(ipar);
1400  }
1401 
1402  // Then release them in the specified order and refit
1403  // --------------------------------------------------
1404  if (debug > 19)
1405  std::cout << " Then release them in order..." << std::endl;
1406 
1407  TString name;
1408  double pval;
1409  double pmin;
1410  double pmax;
1411  double errp;
1412  double errl;
1413  double errh;
1414  int ivar;
1415  double erro;
1416  double cglo;
1417  int n_times = 0;
1418  // n_times = number of loops required to unlock all parameters.
1419 
1420  if (debug > 19)
1421  std::cout << "Before scale parNum" << std::endl;
1422  int scaleParNum = scaleFunction->parNum();
1423  if (debug > 19)
1424  std::cout << "After scale parNum" << std::endl;
1425  // edm::LogInfo("minimizeLikelihood") << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
1426  // edm::LogInfo("minimizeLikelihood") << "number of parameters for resolutionFunction = " << resParNum << std::endl;
1427  // edm::LogInfo("minimizeLikelihood") << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
1428  // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
1429  std::cout << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
1430  std::cout << "number of parameters for resolutionFunction = " << resParNum << std::endl;
1431  std::cout << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
1432  std::cout << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
1433  // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << backgroundFunction->parNum() << std::endl;
1434 
1435  for (int i = 0; i < parnumber; i++) {
1436  // NB ind[] has been set as parorder[] previously
1437  if (n_times < ind[i]) {
1438  edm::LogInfo("minimizeLikelihood") << "n_times = " << n_times << ", ind[" << i << "] = " << ind[i]
1439  << ", scaleParNum = " << scaleParNum << ", doScaleFit[" << loopCounter
1440  << "] = " << doScaleFit[loopCounter] << std::endl;
1441  // Set the n_times only if we will do the fit
1442  if (i < resParNum) {
1443  if (doResolFit[loopCounter])
1444  n_times = ind[i];
1445  } else if (i < resParNum + scaleParNum) {
1446  if (doScaleFit[loopCounter])
1447  n_times = ind[i];
1448  } else if (doBackgroundFit[loopCounter])
1449  n_times = ind[i];
1450  }
1451  }
1452  for (int iorder = 0; iorder < n_times + 1; iorder++) { // Repeat fit n_times times
1453  std::cout << "Starting minimization " << iorder << " of " << n_times << std::endl;
1454 
1455  bool somethingtodo = false;
1456 
1457  // Use parameters from cfg to select which fit to do
1458  // -------------------------------------------------
1459  if (doResolFit[loopCounter]) {
1460  // Release resolution parameters and fit them
1461  // ------------------------------------------
1462  for (unsigned int ipar = 0; ipar < parResol.size(); ++ipar) {
1463  if (parfix[ipar] == 0 && ind[ipar] == iorder) {
1464  rmin.Release(ipar);
1465  somethingtodo = true;
1466  }
1467  }
1468  }
1469  if (doScaleFit[loopCounter]) {
1470  // Release scale parameters and fit them
1471  // -------------------------------------
1472  for (unsigned int ipar = parResol.size(); ipar < parResol.size() + parScale.size(); ++ipar) {
1473  if (parfix[ipar] == 0 && ind[ipar] == iorder) { // parfix=0 means parameter is free
1474  rmin.Release(ipar);
1475  somethingtodo = true;
1476  }
1477  }
1478  scaleFitNotDone_ = false;
1479  }
1480  unsigned int crossSectionParShift = parResol.size() + parScale.size();
1482  // Release cross section parameters and fit them
1483  // ---------------------------------------------
1484  // Note that only cross sections of resonances that are being fitted are released
1485  bool doCrossSection =
1486  crossSectionHandler->releaseParameters(rmin, resfind, parfix, ind, iorder, crossSectionParShift);
1487  if (doCrossSection)
1488  somethingtodo = true;
1489  }
1491  // Release background parameters and fit them
1492  // ------------------------------------------
1493  // for( int ipar=parResol.size()+parScale.size(); ipar<parnumber; ++ipar ) {
1494  // Free only the parameters for the regions, as the resonances intervals are never used to fit the background
1495  unsigned int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
1496  for (unsigned int ipar = bgrParShift; ipar < bgrParShift + backgroundHandler->regionsParNum(); ++ipar) {
1497  // Release only those parameters for the resonances we are fitting
1498  if (parfix[ipar] == 0 && ind[ipar] == iorder &&
1499  backgroundHandler->unlockParameter(resfind, ipar - bgrParShift)) {
1500  rmin.Release(ipar);
1501  somethingtodo = true;
1502  }
1503  }
1504  }
1505 
1506  // OK, now do minimization if some parameter has been released
1507  // -----------------------------------------------------------
1508  if (somethingtodo) {
1509  // #ifdef DEBUG
1510 
1511  std::stringstream fileNum;
1512  fileNum << loopCounter;
1513 
1514  minuitLoop_ = 0;
1515  char name[50];
1516  sprintf(name, "likelihoodInLoop_%d_%d", loopCounter, iorder);
1517  TH1D *tempLikelihoodInLoop = new TH1D(name, "likelihood value in minuit loop", 10000, 0, 10000);
1518  likelihoodInLoop_ = tempLikelihoodInLoop;
1519  char signalProbName[50];
1520  sprintf(signalProbName, "signalProb_%d_%d", loopCounter, iorder);
1521  TH1D *tempSignalProb = new TH1D(signalProbName, "signal probability", 10000, 0, 10000);
1522  signalProb_ = tempSignalProb;
1523  char backgroundProbName[50];
1524  sprintf(backgroundProbName, "backgroundProb_%d_%d", loopCounter, iorder);
1525  TH1D *tempBackgroundProb = new TH1D(backgroundProbName, "background probability", 10000, 0, 10000);
1526  backgroundProb_ = tempBackgroundProb;
1527  // #endif
1528 
1529  // Before we start the minimization we create a list of events with only the events inside a smaller
1530  // window than the one in which the probability is != 0. We will compute the probability for all those
1531  // events and hopefully the margin will avoid them to get a probability = 0 (which causes a discontinuity
1532  // in the likelihood function). The width of this smaller window must be optimized, but we can start using
1533  // an 90% of the normalization window.
1534  double protectionFactor = 0.9;
1535 
1537  for (unsigned int nev = 0; nev < MuScleFitUtils::SavedPair.size(); ++nev) {
1538  const lorentzVector *recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
1539  const lorentzVector *recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
1540  double mass = MuScleFitUtils::invDimuonMass(*recMu1, *recMu2);
1541  // Test all resonances to see if the mass is inside at least one of the windows
1542  bool check = false;
1543  for (int ires = 0; ires < 6; ++ires) {
1544  // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( doBackgroundFit[loopCounter], ires );
1545  std::pair<double, double> windowBorder = backgroundHandler->windowBorders(doBackgroundFit[loopCounter], ires);
1546  // if( resfind[ires] && checkMassWindow( mass, ires, backgroundHandler->resMass( doBackgroundFit[loopCounter], ires ),
1547  // 0.9*windowFactor.first, 0.9*windowFactor.second ) ) {
1548  // double resMassValue = backgroundHandler->resMass( doBackgroundFit[loopCounter], ires );
1549  // double windowBorderLeft = resMassValue - protectionFactor*(resMassValue - windowBorder.first);
1550  // double windowBorderRight = resMassValue + protectionFactor*(windowBorder.second - resMassValue);
1551  double windowBorderShift = (windowBorder.second - windowBorder.first) * (1 - protectionFactor) / 2.;
1552  double windowBorderLeft = windowBorder.first + windowBorderShift;
1553  double windowBorderRight = windowBorder.second - windowBorderShift;
1554  if (resfind[ires] && checkMassWindow(mass, windowBorderLeft, windowBorderRight)) {
1555  check = true;
1556  }
1557  }
1558  if (check) {
1559  MuScleFitUtils::ReducedSavedPair.push_back(std::make_pair(*recMu1, *recMu2));
1560  }
1561  }
1562  std::cout << "Fitting with " << MuScleFitUtils::ReducedSavedPair.size() << " events" << std::endl;
1563 
1564  // rmin.SetMaxIterations(500*parnumber);
1565 
1566  //Print some informations
1567  std::cout << "MINUIT is starting the minimization for the iteration number " << loopCounter << std::endl;
1568 
1569  //Try to set iterations
1570  // rmin.SetMaxIterations(100000);
1571 
1572  std::cout << "maxNumberOfIterations (just set) = " << rmin.GetMaxIterations() << std::endl;
1573 
1575 
1576  // Maximum number of iterations
1577  arglis[0] = 100000;
1578  // tolerance
1579  arglis[1] = 0.1;
1580 
1581  // Run simplex first to get an initial estimate of the minimum
1582  if (startWithSimplex_) {
1583  rmin.mnexcm("SIMPLEX", arglis, 0, ierror);
1584  }
1585 
1586  rmin.mnexcm("MIGRAD", arglis, 2, ierror);
1587 
1588  // #ifdef DEBUG
1589  likelihoodInLoop_->Write();
1590  signalProb_->Write();
1591  backgroundProb_->Write();
1592  delete tempLikelihoodInLoop;
1593  delete tempSignalProb;
1594  delete tempBackgroundProb;
1595  likelihoodInLoop_ = nullptr;
1596  signalProb_ = nullptr;
1597  backgroundProb_ = nullptr;
1598  // #endif
1599 
1600  // Compute again the error matrix
1601  rmin.mnexcm("HESSE", arglis, 0, ierror);
1602 
1603  // Peform minos error analysis.
1604  if (computeMinosErrors_) {
1605  duringMinos_ = true;
1606  rmin.mnexcm("MINOS", arglis, 0, ierror);
1607  duringMinos_ = false;
1608  }
1609 
1610  if (normalizationChanged_ > 1) {
1611  std::cout << "WARNING: normalization changed during fit meaning that events exited from the mass window. This "
1612  "causes a discontinuity in the likelihood function. Please check the scan of the likelihood as a "
1613  "function of the parameters to see if there are discontinuities around the minimum."
1614  << std::endl;
1615  }
1616  }
1617 
1618  // bool notWritten = true;
1619  for (int ipar = 0; ipar < parnumber; ipar++) {
1620  rmin.mnpout(ipar, name, pval, erro, pmin, pmax, ivar);
1621  // Save parameters in parvalue[] vector
1622  // ------------------------------------
1623  if (ierror != 0 && debug > 0) {
1624  std::cout << "[MuScleFitUtils-minimizeLikelihood]: ierror!=0, bogus pars" << std::endl;
1625  }
1626  // for (int ipar=0; ipar<parnumber; ipar++) {
1627  // rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
1628  parvalue[loopCounter][ipar] = pval;
1629  // }
1630 
1631  // int ilax2 = 0;
1632  // Double_t val2pl, val2mi;
1633  // rmin.mnmnot (ipar+1, ilax2, val2pl, val2mi);
1634  rmin.mnerrs(ipar, errh, errl, errp, cglo);
1635 
1636  // Set error on params
1637  // -------------------
1638  if (errp != 0) {
1639  parerr[3 * ipar] = errp;
1640  } else {
1641  parerr[3 * ipar] = (((errh) > (std::abs(errl))) ? (errh) : (std::abs(errl)));
1642  }
1643  parerr[3 * ipar + 1] = errl;
1644  parerr[3 * ipar + 2] = errh;
1645 
1646  if (ipar == 0) {
1647  FitParametersFile << " Resolution fit parameters:" << std::endl;
1648  }
1649  if (ipar == int(parResol.size())) {
1650  FitParametersFile << " Scale fit parameters:" << std::endl;
1651  }
1652  if (ipar == int(parResol.size() + parScale.size())) {
1653  FitParametersFile << " Cross section fit parameters:" << std::endl;
1654  }
1655  if (ipar == int(parResol.size() + parScale.size() + crossSectionHandler->parNum())) {
1656  FitParametersFile << " Background fit parameters:" << std::endl;
1657  }
1658  // if( ipar >= int(parResol.size()+parScale.size()) && ipar < int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) && notWritted ) {
1659 
1660  // std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]));
1661  // std::vector<double>::const_iterator it = relativeCrossSections.begin();
1662  // for( ; it != relativeCrossSections.end(); ++it ) {
1663  // FitParametersFile << " Results of the fit: parameter " << ipar << " has value "
1664  // << *it << "+-" << 0
1665  // << " + " << 0 << " - " << 0
1666  // << " /t/t (" << 0 << ")" << std::endl;
1667  // }
1668 
1669  // notWritten = false;
1670  // }
1671  // else {
1672  FitParametersFile << " Results of the fit: parameter " << ipar << " has value " << pval << "+-"
1673  << parerr[3 * ipar] << " + " << parerr[3 * ipar + 1] << " - "
1674  << parerr[3 * ipar + 2]
1675  // << " \t\t (" << parname[ipar] << ")"
1676  << std::endl;
1677  }
1678  rmin.mnstat(fmin, fdem, errdef, npari, nparx, istat); // NNBB Commented for a check!
1679  FitParametersFile << std::endl;
1680 
1681  if (minimumShapePlots_) {
1682  // Create plots of the minimum vs parameters
1683  // -----------------------------------------
1684  // Keep this after the parameters filling because it recomputes the values and it can compromise the fit results.
1685  if (somethingtodo) {
1686  std::stringstream iorderString;
1687  iorderString << iorder;
1688  std::stringstream iLoopString;
1689  iLoopString << loopCounter;
1690 
1691  for (int ipar = 0; ipar < parnumber; ipar++) {
1692  if (parfix[ipar] == 1)
1693  continue;
1694  std::cout << "plotting parameter = " << ipar + 1 << std::endl;
1695  std::stringstream iparString;
1696  iparString << ipar + 1;
1697  std::stringstream iparStringName;
1698  iparStringName << ipar;
1699  rmin.mncomd(("scan " + iparString.str()).c_str(), ierror);
1700  if (ierror == 0) {
1701  TCanvas *canvas = new TCanvas(("likelihoodCanvas_loop_" + iLoopString.str() + "_oder_" +
1702  iorderString.str() + "_par_" + iparStringName.str())
1703  .c_str(),
1704  ("likelihood_" + iparStringName.str()).c_str(),
1705  1000,
1706  800);
1707  canvas->cd();
1708  // arglis[0] = ipar;
1709  // rmin.mnexcm( "SCA", arglis, 0, ierror );
1710  TGraph *graph = (TGraph *)rmin.GetPlot();
1711  graph->Draw("AP");
1712  // graph->SetTitle(("parvalue["+iparStringName.str()+"]").c_str());
1713  graph->SetTitle(parname[ipar]);
1714  // graph->Write();
1715 
1716  canvas->Write();
1717  }
1718  }
1719 
1720  // // Draw contours of the fit
1721  // TCanvas * canvas = new TCanvas(("contourCanvas_oder_"+iorderString.str()).c_str(), "contour", 1000, 800);
1722  // canvas->cd();
1723  // TGraph * contourGraph = (TGraph*)rmin.Contour(4, 2, 4);
1724  // if( (rmin.GetStatus() == 0) || (rmin.GetStatus() >= 3) ) {
1725  // contourGraph->Draw("AP");
1726  // }
1727  // else {
1728  // std::cout << "Contour graph error: status = " << rmin.GetStatus() << std::endl;
1729  // }
1730  // canvas->Write();
1731  }
1732  }
1733 
1734  } // end loop on iorder
1735  FitParametersFile.close();
1736 
1737  std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters after likelihood " << std::endl;
1738  for (unsigned int ipar = 0; ipar < (unsigned int)parnumber; ipar++) {
1739  std::cout << ipar << " " << parvalue[loopCounter][ipar] << " : free = " << parfix[ipar]
1740  << "; order = " << parorder[ipar] << std::endl;
1741  }
1742 
1743  // Put back parvalue into parResol, parScale, parCrossSection, parBgr
1744  // ------------------------------------------------------------------
1745  for (int i = 0; i < (int)(parResol.size()); ++i) {
1747  }
1748  for (int i = 0; i < (int)(parScale.size()); ++i) {
1749  parScale[i] = parvalue[loopCounter][i + parResol.size()];
1750  }
1751  parCrossSection =
1753  for (unsigned int i = 0; i < parCrossSection.size(); ++i) {
1754  // parCrossSection[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()];
1755  std::cout << "relative cross section[" << i << "] = " << parCrossSection[i] << std::endl;
1756  }
1757  // Save only the fitted background parameters
1758  for (unsigned int i = 0; i < (parBgr.size() - parForResonanceWindows); ++i) {
1760  }
1761 
1762  // Background rescaling from regions to resonances
1763  // -----------------------------------------------
1764  // Only if we fitted the background
1766  // This rescales from regions to resonances
1767  int localMuonType = MuonType;
1768  if (MuonType > 2)
1769  localMuonType = 2;
1771  }
1772 
1773  // Delete the arrays used to set some parameters
1774  delete[] Start;
1775  delete[] Step;
1776  delete[] Mini;
1777  delete[] Maxi;
1778  delete[] ind;
1779  delete[] parname;
1780 }
static std::vector< int > doScaleFit
static std::vector< int > doResolFit
static std::vector< int > parCrossSectionOrder
static std::vector< int > parScaleOrder
static std::vector< int > doCrossSectionFit
int regionsParNum()
Returns the total number of parameters used for the regions.
static std::vector< int > parBgrOrder
static std::vector< int > parfix
static unsigned int loopCounter
static std::vector< double > parResolMax
void likelihood(int &npar, double *grad, double &fval, double *xval, int flag)
static bool startWithSimplex_
static std::vector< int > doBackgroundFit
static std::vector< double > parResol
void setParameters(double *Start, double *Step, double *Mini, double *Maxi, int *ind, TString *parname, const std::vector< double > &parCrossSection, const std::vector< int > &parCrossSectionOrder, const std::vector< int > &resfind)
Initializes the arrays needed by Minuit.
static TH1D * backgroundProb_
static BackgroundHandler * backgroundHandler
static int debug
static double ResMass[6]
static double massWindowHalfWidth[3][6]
static bool scaleFitNotDone_
static unsigned int normalizationChanged_
static bool checkMassWindow(const double &mass, const double &leftBorder, const double &rightBorder)
Method to check if the mass value is within the mass window of the i-th resonance.
static std::vector< int > parBgrFix
static std::vector< double > parResolMin
static bool minimumShapePlots_
static std::vector< int > parScaleFix
static bool computeMinosErrors_
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
static scaleFunctionBase< std::vector< double > > * scaleFunctionForVec
static std::vector< std::vector< double > > parvalue
static int ScaleFitType
static std::vector< std::pair< lorentzVector, lorentzVector > > ReducedSavedPair
static int minuitLoop_
static std::vector< double > parScaleMin
void rescale(std::vector< double > &parBgr, const double *ResMass, const double *massWindowHalfWidth, const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &muonPairs, const double &weight=1.)
static std::vector< double > parBgr
virtual void setParameters(double *Start, double *Step, double *Mini, double *Maxi, int *ind, TString *parname, const T &parScale, const std::vector< int > &parScaleOrder, const int muonType)=0
This method is used to differentiate parameters among the different functions.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool releaseParameters(TMinuit &rmin, const std::vector< int > &resfind, const std::vector< int > &parfix, const int *ind, const int iorder, const unsigned int shift)
Use the information in resfind, parorder and parfix to release the N-1 variables. ...
static std::vector< double > parScaleStep
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
static double invDimuonMass(const lorentzVector &mu1, const lorentzVector &mu2)
static TH1D * signalProb_
static std::vector< int > parResolFix
static TMinuit * rminPtr_
static std::vector< int > parResolOrder
static int BgrFitType
virtual int parNum() const
Definition: Functions.h:76
static TH1D * likelihoodInLoop_
Log< level::Info, false > LogInfo
static std::vector< double > parScaleMax
static std::vector< double > parScale
virtual int parNum() const
Definition: Functions.h:865
std::vector< double > relativeCrossSections(const double *variables, const std::vector< int > &resfind)
Perform a variable transformation from N-1 to relative cross sections.
void setParameters(double *Start, double *Step, double *Mini, double *Maxi, int *ind, TString *parname, const std::vector< double > &parBgr, const std::vector< int > &parBgrOrder, const int muonType)
Sets initial parameters for all the functions.
virtual void resetParameters(std::vector< double > *scaleVec) const
This method is used to reset the scale parameters to neutral values (useful for iterations > 0) ...
Definition: Functions.h:46
static bool duringMinos_
static int ResolFitType
static std::vector< double > parCrossSection
std::pair< double, double > windowBorders(const bool doBackgroundFit, const int ires)
Returns the appropriate window borders depending on whether the background is being fitted and on the...
bool unlockParameter(const std::vector< int > &resfind, const unsigned int ipar)
returns true if the parameter is to be unlocked
def canvas(sub, attr)
Definition: svgfig.py:482
virtual void setParameters(double *Start, double *Step, double *Mini, double *Maxi, int *ind, TString *parname, const T &parResol, const std::vector< int > &parResolOrder, const int muonType)
This method is used to differentiate parameters among the different functions.
Definition: Functions.h:841
static int MuonType
static scaleFunctionBase< double * > * scaleFunction
static std::vector< int > resfind
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
static std::vector< int > parorder
static int FitStrategy
static CrossSectionHandler * crossSectionHandler
static std::vector< double > parResolStep

◆ probability()

double MuScleFitUtils::probability ( const double &  mass,
const double &  massResol,
const double  GLvalue[][1001][1001],
const double  GLnorm[][1001],
const int  iRes,
const int  iY 
)
static

Computes the probability given the mass, mass resolution and the arrays with the probabilities and the normalizations.

After the introduction of the rapidity bins for the Z the probability method works in the following way:

  • if passing iRes == 0, iY is used to select the rapidity bin
  • if passing iRes != 0, iY is used to select the resonance

Definition at line 750 of file MuScleFitUtils.cc.

References counter_resprob, gather_cfg::cout, debug, createfilelist::int, LogDebug, EgHLTOffHistBins_cfi::mass, nbins, ResHalfWidth, ResMass, ResMaxSigma, and ResMinMass.

Referenced by massProb(), and MuScleFitBase::ProbForIntegral::operator()().

755  {
756  if (iRes == 0 && iY > 23) {
757  std::cout << "WARNING: rapidity bin selected = " << iY << " but there are only histograms for the first 24 bins"
758  << std::endl;
759  }
760 
761  double PS = 0.;
762  bool insideProbMassWindow = true;
763  // Interpolate the four values of GLZValue[] in the
764  // grid square within which the (mass,sigma) values lay
765  // ----------------------------------------------------
766  // This must be done with respect to the width used in the computation of the probability distribution,
767  // so that the bin 0 really matches the bin 0 of that distribution.
768  // double fracMass = (mass-(ResMass[iRes]-ResHalfWidth[iRes]))/(2*ResHalfWidth[iRes]);
769  double fracMass = (mass - ResMinMass[iRes]) / (2 * ResHalfWidth[iRes]);
770  if (debug > 1)
771  std::cout << std::setprecision(9) << "mass ResMinMass[iRes] ResHalfWidth[iRes] ResHalfWidth[iRes]" << mass << " "
772  << ResMinMass[iRes] << " " << ResHalfWidth[iRes] << " " << ResHalfWidth[iRes] << std::endl;
773  int iMassLeft = (int)(fracMass * (double)nbins);
774  int iMassRight = iMassLeft + 1;
775  double fracMassStep = (double)nbins * (fracMass - (double)iMassLeft / (double)nbins);
776  if (debug > 1)
777  std::cout << "nbins iMassLeft fracMass " << nbins << " " << iMassLeft << " " << fracMass << std::endl;
778 
779  // Simple protections for the time being: the region where we fit should not include
780  // values outside the boundaries set by ResMass-ResHalfWidth : ResMass+ResHalfWidth
781  // ---------------------------------------------------------------------------------
782  if (iMassLeft < 0) {
783  edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassLeft=" << iMassLeft
784  << "; mass = " << mass << " and bounds are " << ResMinMass[iRes] << ":"
785  << ResMinMass[iRes] + 2 * ResHalfWidth[iRes] << " - iMassLeft set to 0" << std::endl;
786  iMassLeft = 0;
787  iMassRight = 1;
788  insideProbMassWindow = false;
789  }
790  if (iMassRight > nbins) {
791  edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassRight=" << iMassRight
792  << "; mass = " << mass << " and bounds are " << ResMinMass[iRes] << ":"
793  << ResMass[iRes] + 2 * ResHalfWidth[iRes] << " - iMassRight set to " << nbins - 1
794  << std::endl;
795  iMassLeft = nbins - 1;
796  iMassRight = nbins;
797  insideProbMassWindow = false;
798  }
799  double fracSigma = (massResol / ResMaxSigma[iRes]);
800  int iSigmaLeft = (int)(fracSigma * (double)nbins);
801  int iSigmaRight = iSigmaLeft + 1;
802  double fracSigmaStep = (double)nbins * (fracSigma - (double)iSigmaLeft / (double)nbins);
803  // std::cout << "massResol = " << massResol << std::endl;
804  // std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
805  // std::cout << "fracSigma = " << fracSigma << std::endl;
806  // std::cout << "nbins = " << nbins << std::endl;
807  // std::cout << "ISIGMALEFT = " << iSigmaLeft << std::endl;
808  // std::cout << "ISIGMARIGHT = " << iSigmaRight << std::endl;
809  // std::cout << "fracSigmaStep = " << fracSigmaStep << std::endl;
810 
811  // Simple protections for the time being: they should not affect convergence, since
812  // ResMaxSigma is set to very large values, and if massResol exceeds them the fit
813  // should not get any prize for that (for large sigma, the prob. distr. becomes flat)
814  // ----------------------------------------------------------------------------------
815  if (iSigmaLeft < 0) {
816  edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaLeft=" << iSigmaLeft
817  << ", with massResol = " << massResol
818  << " and ResMaxSigma[iRes] = " << ResMaxSigma[iRes] << " - iSigmaLeft set to 0"
819  << std::endl;
820  iSigmaLeft = 0;
821  iSigmaRight = 1;
822  }
823  if (iSigmaRight > nbins) {
824  if (counter_resprob < 100)
825  edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaRight=" << iSigmaRight
826  << ", with massResol = " << massResol
827  << " and ResMaxSigma[iRes] = " << ResMaxSigma[iRes] << " - iSigmaRight set to "
828  << nbins - 1 << std::endl;
829  iSigmaLeft = nbins - 1;
830  iSigmaRight = nbins;
831  }
832 
833  // If f11,f12,f21,f22 are the values at the four corners, one finds by linear interpolation the
834  // formula below for PS
835  // --------------------------------------------------------------------------------------------
836  if (insideProbMassWindow) {
837  double f11 = 0.;
838  if (GLnorm[iY][iSigmaLeft] > 0)
839  f11 = GLvalue[iY][iMassLeft][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
840  double f12 = 0.;
841  if (GLnorm[iY][iSigmaRight] > 0)
842  f12 = GLvalue[iY][iMassLeft][iSigmaRight] / GLnorm[iY][iSigmaRight];
843  double f21 = 0.;
844  if (GLnorm[iY][iSigmaLeft] > 0)
845  f21 = GLvalue[iY][iMassRight][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
846  double f22 = 0.;
847  if (GLnorm[iY][iSigmaRight] > 0)
848  f22 = GLvalue[iY][iMassRight][iSigmaRight] / GLnorm[iY][iSigmaRight];
849  PS = f11 + (f12 - f11) * fracSigmaStep + (f21 - f11) * fracMassStep +
850  (f22 - f21 - f12 + f11) * fracMassStep * fracSigmaStep;
851  if (PS > 0.1 || debug > 1)
852  LogDebug("MuScleFitUtils") << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22=" << f11 << " " << f12 << " "
853  << f21 << " " << f22 << " "
854  << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR=" << iSigmaLeft
855  << " " << iSigmaRight << " GLvalue[" << iY << "][" << iMassLeft
856  << "] = " << GLvalue[iY][iMassLeft][iSigmaLeft] << " GLnorm[" << iY << "]["
857  << iSigmaLeft << "] = " << GLnorm[iY][iSigmaLeft] << std::endl;
858 
859  // if (PS>0.1) std::cout << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
860  // << f11 << " " << f12 << " " << f21 << " " << f22 << " "
861  // << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
862  // << iSigmaLeft << " " << iSigmaRight << " GLV,GLN="
863  // << GLvalue[iY][iMassLeft][iSigmaLeft]
864  // << " " << GLnorm[iY][iSigmaLeft] << std::endl;
865 
866  } else {
867  edm::LogInfo("probability") << "outside mass probability window. Setting PS[" << iRes << "] = 0" << std::endl;
868  }
869 
870  // if( PS != PS ) {
871  // std::cout << "ERROR: PS = " << PS << " for iRes = " << iRes << std::endl;
872 
873  // std::cout << "mass = " << mass << ", massResol = " << massResol << std::endl;
874  // std::cout << "fracMass = " << fracMass << ", iMassLeft = " << iMassLeft
875  // << ", iMassRight = " << iMassRight << ", fracMassStep = " << fracMassStep << std::endl;
876  // std::cout << "fracSigma = " << fracSigma << ", iSigmaLeft = " << iSigmaLeft
877  // << ", iSigmaRight = " << iSigmaRight << ", fracSigmaStep = " << fracSigmaStep << std::endl;
878  // std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
879 
880  // std::cout << "GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
881  // << " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;
882  // }
883 
884  return PS;
885 }
static int nbins
static double ResMinMass[6]
static int debug
static double ResMass[6]
static double ResMaxSigma[6]
Log< level::Info, false > LogInfo
static double ResHalfWidth[6]
static int counter_resprob
#define LogDebug(id)

Member Data Documentation

◆ backgroundFunctionsRegions

const int MuScleFitUtils::backgroundFunctionsRegions
static

Definition at line 163 of file MuScleFitUtils.h.

◆ backgroundHandler

BackgroundHandler * MuScleFitUtils::backgroundHandler
static

◆ backgroundProb_

TH1D * MuScleFitUtils::backgroundProb_ = nullptr
static

Definition at line 186 of file MuScleFitUtils.h.

Referenced by massProb(), and minimizeLikelihood().

◆ BgrFitType

int MuScleFitUtils::BgrFitType = 0
static

Definition at line 158 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood().

◆ biasFunction

scaleFunctionBase< std::vector< double > > * MuScleFitUtils::biasFunction = nullptr
static

Definition at line 151 of file MuScleFitUtils.h.

Referenced by applyBias(), and MuScleFit::MuScleFit().

◆ BiasType

int MuScleFitUtils::BiasType = 0
static

◆ computeMinosErrors_

bool MuScleFitUtils::computeMinosErrors_
static

Definition at line 280 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ counter_resprob

int MuScleFitUtils::counter_resprob = 0
static

Definition at line 215 of file MuScleFitUtils.h.

Referenced by massResolution(), probability(), and MuScleFit::startingNewLoop().

◆ crossSection

double MuScleFitUtils::crossSection[6]
static

Definition at line 138 of file MuScleFitUtils.h.

◆ crossSectionHandler

CrossSectionHandler * MuScleFitUtils::crossSectionHandler
static

◆ debug

int MuScleFitUtils::debug = 0
static

◆ debugMassResol_

bool MuScleFitUtils::debugMassResol_
static

◆ deltaPhiMaxCut_

double MuScleFitUtils::deltaPhiMaxCut_ = 100.
static

Definition at line 266 of file MuScleFitUtils.h.

Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ deltaPhiMinCut_

double MuScleFitUtils::deltaPhiMinCut_ = -100.
static

Definition at line 265 of file MuScleFitUtils.h.

Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ doBackgroundFit

std::vector< int > MuScleFitUtils::doBackgroundFit
static

◆ doCrossSectionFit

std::vector< int > MuScleFitUtils::doCrossSectionFit
static

◆ doResolFit

std::vector< int > MuScleFitUtils::doResolFit
static

◆ doScaleFit

std::vector< int > MuScleFitUtils::doScaleFit
static

◆ duringMinos_

bool MuScleFitUtils::duringMinos_ = false
static

Definition at line 188 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood().

◆ FitStrategy

int MuScleFitUtils::FitStrategy = 1
static

Definition at line 211 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ genMuscleFitPair

std::vector< std::pair< MuScleFitMuon, MuScleFitMuon > > MuScleFitUtils::genMuscleFitPair
static

Definition at line 235 of file MuScleFitUtils.h.

Referenced by MuScleFit::selectMuons().

◆ genPair

std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::genPair
static

Definition at line 233 of file MuScleFitUtils.h.

Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().

◆ GLNorm

double MuScleFitUtils::GLNorm
static

◆ GLValue

double MuScleFitUtils::GLValue
static

◆ GLZNorm

double MuScleFitUtils::GLZNorm
static

◆ GLZValue

double MuScleFitUtils::GLZValue
static

◆ goodmuon

int MuScleFitUtils::goodmuon = 0
static

◆ iev_

int MuScleFitUtils::iev_ = 0
static

◆ likelihoodInLoop_

TH1D * MuScleFitUtils::likelihoodInLoop_ = nullptr
static

Definition at line 184 of file MuScleFitUtils.h.

Referenced by likelihood(), and minimizeLikelihood().

◆ loopCounter

unsigned int MuScleFitUtils::loopCounter = 5
static

◆ massResolComponents

MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents
static

Definition at line 262 of file MuScleFitUtils.cc.

Referenced by MuScleFit::duringFastLoop(), and massResolution().

◆ massWindowHalfWidth

double MuScleFitUtils::massWindowHalfWidth
static

Definition at line 134 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ maxMuonEtaFirstRange_

double MuScleFitUtils::maxMuonEtaFirstRange_ = 6.
static

Definition at line 262 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ maxMuonEtaSecondRange_

double MuScleFitUtils::maxMuonEtaSecondRange_ = 100.
static

Definition at line 264 of file MuScleFitUtils.h.

Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ maxMuonPt_

double MuScleFitUtils::maxMuonPt_ = 100000000.
static

Definition at line 260 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ minimumShapePlots_

bool MuScleFitUtils::minimumShapePlots_
static

Definition at line 281 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ minMuonEtaFirstRange_

double MuScleFitUtils::minMuonEtaFirstRange_ = -6.
static

Definition at line 261 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ minMuonEtaSecondRange_

double MuScleFitUtils::minMuonEtaSecondRange_ = -100.
static

Definition at line 263 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ minMuonPt_

double MuScleFitUtils::minMuonPt_ = 0.
static

Definition at line 259 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ minuitLoop_

int MuScleFitUtils::minuitLoop_ = 0
static

Definition at line 183 of file MuScleFitUtils.h.

Referenced by likelihood(), massProb(), and minimizeLikelihood().

◆ mMu2

const double MuScleFitUtils::mMu2 = 0.011163612
static

◆ motherPdgIdArray

const unsigned int MuScleFitUtils::motherPdgIdArray = {23, 100553, 100553, 553, 100443, 443}
static

Definition at line 143 of file MuScleFitUtils.h.

Referenced by findGenMuFromRes(), and findSimMuFromRes().

◆ muMass

const double MuScleFitUtils::muMass = 0.105658
static

Definition at line 140 of file MuScleFitUtils.h.

Referenced by fromPtEtaPhiToPxPyPz().

◆ MuonType

int MuScleFitUtils::MuonType
static

Definition at line 223 of file MuScleFitUtils.h.

Referenced by massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ MuonTypeForCheckMassWindow

int MuScleFitUtils::MuonTypeForCheckMassWindow
static

Definition at line 224 of file MuScleFitUtils.h.

Referenced by MuScleFit::MuScleFit().

◆ nbins

int MuScleFitUtils::nbins = 1000
static

◆ normalizationChanged_

unsigned int MuScleFitUtils::normalizationChanged_ = 0
static

Definition at line 245 of file MuScleFitUtils.h.

Referenced by likelihood(), and minimizeLikelihood().

◆ normalizeLikelihoodByEventNumber_

bool MuScleFitUtils::normalizeLikelihoodByEventNumber_ = true
static

Definition at line 240 of file MuScleFitUtils.h.

Referenced by likelihood(), and MuScleFit::MuScleFit().

◆ oldNormalization_

double MuScleFitUtils::oldNormalization_ = 0.
static

Definition at line 244 of file MuScleFitUtils.h.

Referenced by likelihood(), and MuScleFit::startingNewLoop().

◆ parBgr

std::vector< double > MuScleFitUtils::parBgr
static

◆ parBgrFix

std::vector< int > MuScleFitUtils::parBgrFix
static

◆ parBgrOrder

std::vector< int > MuScleFitUtils::parBgrOrder
static

◆ parBias

std::vector< double > MuScleFitUtils::parBias
static

Definition at line 191 of file MuScleFitUtils.h.

Referenced by applyBias(), MuScleFit::checkParameters(), and MuScleFit::MuScleFit().

◆ parCrossSection

std::vector< double > MuScleFitUtils::parCrossSection
static

◆ parCrossSectionFix

std::vector< int > MuScleFitUtils::parCrossSectionFix
static

Definition at line 204 of file MuScleFitUtils.h.

Referenced by MuScleFit::checkParameters(), and MuScleFit::MuScleFit().

◆ parCrossSectionOrder

std::vector< int > MuScleFitUtils::parCrossSectionOrder
static

◆ parfix

std::vector<int> MuScleFitUtils::parfix
static

Definition at line 228 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood().

◆ parorder

std::vector<int> MuScleFitUtils::parorder
static

Definition at line 229 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood().

◆ parResol

std::vector< double > MuScleFitUtils::parResol
static

◆ parResolFix

std::vector< int > MuScleFitUtils::parResolFix
static

◆ parResolMax

std::vector< double > MuScleFitUtils::parResolMax
static

Definition at line 195 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parResolMin

std::vector< double > MuScleFitUtils::parResolMin
static

Definition at line 194 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parResolOrder

std::vector< int > MuScleFitUtils::parResolOrder
static

◆ parResolStep

std::vector< double > MuScleFitUtils::parResolStep
static

Definition at line 193 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parScale

std::vector< double > MuScleFitUtils::parScale
static

◆ parScaleFix

std::vector< int > MuScleFitUtils::parScaleFix
static

◆ parScaleMax

std::vector< double > MuScleFitUtils::parScaleMax
static

Definition at line 199 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parScaleMin

std::vector< double > MuScleFitUtils::parScaleMin
static

Definition at line 198 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parScaleOrder

std::vector< int > MuScleFitUtils::parScaleOrder
static

◆ parScaleStep

std::vector< double > MuScleFitUtils::parScaleStep
static

Definition at line 197 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ parSmear

std::vector< double > MuScleFitUtils::parSmear
static

◆ parvalue

std::vector< std::vector< double > > MuScleFitUtils::parvalue
static

Definition at line 226 of file MuScleFitUtils.h.

Referenced by MuScleFit::duringFastLoop(), and minimizeLikelihood().

◆ rapidityBinsForZ_

bool MuScleFitUtils::rapidityBinsForZ_ = true
static

◆ ReducedSavedPair

std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::ReducedSavedPair
static

Definition at line 232 of file MuScleFitUtils.h.

Referenced by likelihood(), and minimizeLikelihood().

◆ resfind

std::vector< int > MuScleFitUtils::resfind
static

◆ ResFound

bool MuScleFitUtils::ResFound = false
static

◆ ResGamma

double MuScleFitUtils::ResGamma = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932}
static

Definition at line 135 of file MuScleFitUtils.h.

Referenced by massProb().

◆ ResHalfWidth

double MuScleFitUtils::ResHalfWidth
static

◆ ResMass

double MuScleFitUtils::ResMass = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969}
static

◆ ResMaxSigma

double MuScleFitUtils::ResMaxSigma
static

◆ ResMinMass

double MuScleFitUtils::ResMinMass = {-99, -99, -99, -99, -99, -99}
static

◆ ResolFitType

int MuScleFitUtils::ResolFitType = 0
static

◆ resolutionFunction

resolutionFunctionBase< double * > * MuScleFitUtils::resolutionFunction = nullptr
static

◆ resolutionFunctionForVec

resolutionFunctionBase< std::vector< double > > * MuScleFitUtils::resolutionFunctionForVec = nullptr
static

◆ rminPtr_

TMinuit * MuScleFitUtils::rminPtr_ = nullptr
static

Definition at line 242 of file MuScleFitUtils.h.

Referenced by likelihood(), and minimizeLikelihood().

◆ SavedPair

std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::SavedPair
static

◆ SavedPairMuScleFitMuons

std::vector< std::pair< MuScleFitMuon, MuScleFitMuon > > MuScleFitUtils::SavedPairMuScleFitMuons
static

Definition at line 234 of file MuScleFitUtils.h.

Referenced by MuScleFit::selectMuons().

◆ scaleFitNotDone_

bool MuScleFitUtils::scaleFitNotDone_ = true
static

Definition at line 238 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood().

◆ ScaleFitType

int MuScleFitUtils::ScaleFitType = 0
static

Definition at line 155 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ scaleFunction

scaleFunctionBase< double * > * MuScleFitUtils::scaleFunction = nullptr
static

Definition at line 156 of file MuScleFitUtils.h.

Referenced by applyScale(), minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ scaleFunctionForVec

scaleFunctionBase< std::vector< double > > * MuScleFitUtils::scaleFunctionForVec = nullptr
static

Definition at line 157 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ separateRanges_

bool MuScleFitUtils::separateRanges_ = true
static

Definition at line 258 of file MuScleFitUtils.h.

Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().

◆ sherpa_

bool MuScleFitUtils::sherpa_ = false
static

Definition at line 248 of file MuScleFitUtils.h.

Referenced by findGenMuFromRes(), and MuScleFit::MuScleFit().

◆ signalProb_

TH1D * MuScleFitUtils::signalProb_ = nullptr
static

Definition at line 185 of file MuScleFitUtils.h.

Referenced by massProb(), and minimizeLikelihood().

◆ simPair

std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::simPair
static

Definition at line 236 of file MuScleFitUtils.h.

Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().

◆ smearFunction

smearFunctionBase * MuScleFitUtils::smearFunction = nullptr
static

Definition at line 148 of file MuScleFitUtils.h.

Referenced by applySmearing(), and MuScleFit::MuScleFit().

◆ SmearType

int MuScleFitUtils::SmearType = 0
static

◆ speedup

bool MuScleFitUtils::speedup = false
static

◆ startWithSimplex_

bool MuScleFitUtils::startWithSimplex_
static

Definition at line 279 of file MuScleFitUtils.h.

Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().

◆ totalResNum

const int MuScleFitUtils::totalResNum = 6
static

Definition at line 133 of file MuScleFitUtils.h.

Referenced by massProb().

◆ useProbsFile_

bool MuScleFitUtils::useProbsFile_ = true
static

◆ x

double MuScleFitUtils::x
static