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, MuScleFitMuon > findBestRecoRes (const std::vector< MuScleFitMuon > &muons)
 
static std::pair< SimTrack, SimTrack > findBestSimuRes (const std::vector< SimTrack > &simMuons)
 
static std::pair< lorentzVector, lorentzVector > findGenMuFromRes (const edm::HepMCProduct *evtMC)
 
static std::pair< lorentzVector, lorentzVector > findGenMuFromRes (const reco::GenParticleCollection *genParticles)
 
static std::pair< lorentzVector, lorentzVector > findSimMuFromRes (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 ResolutionFunction &resolFunc)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, double *parval)
 
static double massResolution (const lorentzVector &mu1, const lorentzVector &mu2, std::unique_ptr< double > parval)
 
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 BackgroundHandler * backgroundHandler
 
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 CrossSectionHandler * crossSectionHandler
 
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 smearFunctionBase * smearFunction = 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.

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 }

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

Referenced by MuScleFit::applyBias().

◆ 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.

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 }

References chg, triggerObjects_cff::id, createfilelist::int, and AlCaHLTBitMon_ParallelJobs::p.

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

◆ applyScale() [2/2]

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

Definition at line 483 of file MuScleFitUtils.cc.

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 }

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

◆ applySmearing()

lorentzVector MuScleFitUtils::applySmearing ( const lorentzVector &  muon)
static

Definition at line 416 of file MuScleFitUtils.cc.

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
428  smearFunction->smear(pt, eta, phi, y, parSmear);
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 }

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

◆ 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.

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

References EgHLTOffHistBins_cfi::mass.

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

◆ computeWeight()

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

Definition at line 1163 of file MuScleFitUtils.cc.

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 }

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

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

◆ deltaPhi()

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

Definition at line 109 of file MuScleFitUtils.h.

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  }

References Pi.

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

◆ 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.

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  }

References deltaPhi(), and Pi.

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

◆ 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.

126  {
127  return sqrt(std::pow(eta1 - eta2, 2) + std::pow(deltaPhi(phi1, phi2), 2));
128  }

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

◆ findBestRecoRes()

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

Definition at line 315 of file MuScleFitUtils.cc.

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 }

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

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

◆ findBestSimuRes()

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

Definition at line 279 of file MuScleFitUtils.cc.

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 }

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

Referenced by MuScleFitPlotter::fillSim().

◆ findGenMuFromRes() [1/2]

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

Definition at line 2306 of file MuScleFitUtils.cc.

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 }

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

◆ findGenMuFromRes() [2/2]

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

Definition at line 2344 of file MuScleFitUtils.cc.

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 }

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

Referenced by MuScleFitGenFilter::filter().

◆ 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.

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 }

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

Referenced by MuScleFitPlotter::fillGenSim().

◆ fitMass()

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

Definition at line 1985 of file MuScleFitUtils.cc.

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 }

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

◆ fitReso()

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

Definition at line 2098 of file MuScleFitUtils.cc.

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 }

References hltPixelTracks_cff::chi2, generateEDF::cont, gather_cfg::cout, MillePedeFileConverter_cfg::e, submitPVResolutionJobs::err, Gaussian(), timingPdfMaker::histo, mps_fire::i, dqmiolumiharvest::j, Skims_PA_cff::name, groupFilesInBlocks::nn, bookConverter::results, x, and geometryCSVtoXML::xx.

◆ fromPtEtaPhiToPxPyPz()

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

Definition at line 502 of file MuScleFitUtils.cc.

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 }

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

◆ invDimuonMass()

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

Definition at line 518 of file MuScleFitUtils.cc.

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

Referenced by likelihood(), and minimizeLikelihood().

◆ 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.

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 }

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

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

◆ 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.

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 }

References HLT_FULL_cff::eta1, HLT_FULL_cff::eta2, triggerObjects_cff::id, createfilelist::int, EgHLTOffHistBins_cfi::mass, massProb(), and AlCaHLTBitMon_ParallelJobs::p.

◆ 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.

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
1000  if (MuScleFitUtils::rapidityBinsForZ_) {
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]),
1032  MuScleFitUtils::totalResNum,
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;
1070  if (!MuScleFitUtils::rapidityBinsForZ_)
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 =
1085  backgroundHandler->backgroundFunction(doBackgroundFit[loopCounter],
1086  &(parval[bgrParShift]),
1087  MuScleFitUtils::totalResNum,
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(
1138  MuScleFitUtils::minuitLoop_,
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(
1146  MuScleFitUtils::minuitLoop_, MuScleFitUtils::backgroundProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PB);
1147  }
1148  return P;
1149 }

References funct::abs(), BackgroundHandler::backgroundFunction(), backgroundHandler, backgroundProb_, checkMassWindow(), gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, HLT_FULL_cff::eta1, HLT_FULL_cff::eta2, GLNorm, GLValue, GLZNorm, GLZValue, mps_fire::i, createfilelist::int, ires, 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().

◆ massResolution() [1/5]

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

Referenced by ResolutionAnalyzer::analyze(), MuScleFit::duringFastLoop(), ErrorsAnalyzer::fillHistograms(), findBestRecoRes(), likelihood(), and massResolution().

◆ massResolution() [2/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.

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 }

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

◆ massResolution() [3/5]

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

Definition at line 524 of file MuScleFitUtils.cc.

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 }

References triggerObjects_cff::id, createfilelist::int, massResolution(), and AlCaHLTBitMon_ParallelJobs::p.

◆ 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.

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 }

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_FULL_cff::eta1, HLT_FULL_cff::eta2, JetChargeProducer_cfi::exp, ires, LogDebug, EgHLTOffHistBins_cfi::mass, massResolComponents, mMu2, funct::pow(), HLT_FULL_cff::pt1, HLT_FULL_cff::pt2, resfind, ResHalfWidth, ResMass, ResMaxSigma, resolutionFunction, resolutionFunctionBase< T >::sigmaCotgTh(), resolutionFunctionBase< T >::sigmaPhi(), resolutionFunctionBase< T >::sigmaPt(), funct::sin(), and mathSSE::sqrt().

◆ massResolution() [5/5]

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

◆ minimizeLikelihood()

void MuScleFitUtils::minimizeLikelihood ( )
static

Definition at line 1197 of file MuScleFitUtils.cc.

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()) {
1311  MuScleFitUtils::resolutionFunctionForVec->setParameters(
1312  Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parResolStep, parResolMin, parResolMax, MuonType);
1313  } else {
1314  MuScleFitUtils::resolutionFunctionForVec->setParameters(
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
1319  int resParNum = MuScleFitUtils::resolutionFunctionForVec->parNum();
1320 
1321  if (!parScaleStep.empty() && !parScaleMin.empty() && !parScaleMax.empty()) {
1322  MuScleFitUtils::scaleFunctionForVec->setParameters(&(Start[resParNum]),
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 {
1335  MuScleFitUtils::scaleFunctionForVec->setParameters(&(Start[resParNum]),
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]),
1354  parCrossSection,
1355  parCrossSectionOrder,
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();
1481  if (doCrossSectionFit[loopCounter]) {
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  }
1490  if (doBackgroundFit[loopCounter]) {
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 
1536  MuScleFitUtils::ReducedSavedPair.clear();
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 
1574  MuScleFitUtils::normalizationChanged_ = 0;
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) {
1746  parResol[i] = parvalue[loopCounter][i];
1747  }
1748  for (int i = 0; i < (int)(parScale.size()); ++i) {
1749  parScale[i] = parvalue[loopCounter][i + parResol.size()];
1750  }
1751  parCrossSection =
1752  crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size() + parScale.size()]), resfind);
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) {
1759  parBgr[i] = parvalue[loopCounter][i + parResol.size() + parScale.size() + crossSectionHandler->parNum()];
1760  }
1761 
1762  // Background rescaling from regions to resonances
1763  // -----------------------------------------------
1764  // Only if we fitted the background
1765  if (doBackgroundFit[loopCounter]) {
1766  // This rescales from regions to resonances
1767  int localMuonType = MuonType;
1768  if (MuonType > 2)
1769  localMuonType = 2;
1770  backgroundHandler->rescale(parBgr, ResMass, massWindowHalfWidth[localMuonType], MuScleFitUtils::ReducedSavedPair);
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 }

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(), iorder, ires, likelihood(), likelihoodInLoop_, loopCounter, EgHLTOffHistBins_cfi::mass, massWindowHalfWidth, minimumShapePlots_, minuitLoop_, MuonType, Skims_PA_cff::name, 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, trackingParticleMuon_cfi::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(), CrossSectionHandler::setParameters(), BackgroundHandler::setParameters(), resolutionFunctionBase< T >::setParameters(), signalProb_, startWithSimplex_, BackgroundHandler::unlockParameter(), and BackgroundHandler::windowBorders().

Referenced by MuScleFit::endOfFastLoop().

◆ 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.

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 }

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

Member Data Documentation

◆ backgroundFunctionsRegions

const int MuScleFitUtils::backgroundFunctionsRegions
static

Definition at line 163 of file MuScleFitUtils.h.

◆ backgroundHandler

BackgroundHandler * MuScleFitUtils::backgroundHandler
static

Definition at line 175 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 149 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 171 of file MuScleFitUtils.h.

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

◆ debug

int MuScleFitUtils::debug = 0
static

Definition at line 130 of file MuScleFitUtils.h.

Referenced by applyBias(), applyScale(), applySmearing(), computeWeight(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), findBestRecoRes(), findGenMuFromRes(), fitMass(), likelihood(), runTauIdMVA.TauIDEmbedder::loadMVA_WPs_run2_2017(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), probability(), and runTauIdMVA.TauIDEmbedder::runTauID().

◆ debugMassResol_

bool MuScleFitUtils::debugMassResol_
static

Definition at line 268 of file MuScleFitUtils.h.

Referenced by MuScleFit::duringFastLoop(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), MuScleFitBase::fillHistoMap(), massResolution(), and MuScleFit::MuScleFit().

◆ 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

Definition at line 181 of file MuScleFitUtils.h.

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

◆ doCrossSectionFit

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

Definition at line 180 of file MuScleFitUtils.h.

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

◆ doResolFit

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

Definition at line 178 of file MuScleFitUtils.h.

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

◆ doScaleFit

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

Definition at line 179 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 219 of file MuScleFitUtils.h.

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

◆ GLValue

double MuScleFitUtils::GLValue
static

Definition at line 218 of file MuScleFitUtils.h.

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

◆ GLZNorm

double MuScleFitUtils::GLZNorm
static

Definition at line 217 of file MuScleFitUtils.h.

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

◆ GLZValue

double MuScleFitUtils::GLZValue
static

Definition at line 216 of file MuScleFitUtils.h.

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

◆ goodmuon

int MuScleFitUtils::goodmuon = 0
static

Definition at line 214 of file MuScleFitUtils.h.

Referenced by MuScleFit::applyBias(), applySmearing(), MuScleFit::applySmearing(), ResolutionAnalyzer::fillMuonCollection(), MuScleFit::fillMuonCollection(), MuScleFit::MuScleFit(), and MuScleFit::startingNewLoop().

◆ iev_

int MuScleFitUtils::iev_ = 0
static

Definition at line 253 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 145 of file MuScleFitUtils.h.

Referenced by computeWeight(), likelihood(), massProb(), minimizeLikelihood(), and MuScleFit::startingNewLoop().

◆ 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

Definition at line 139 of file MuScleFitUtils.h.

Referenced by ResolutionAnalyzer::analyze(), ResolutionAnalyzer::fillMuonCollection(), MuScleFit::fillMuonCollection(), ErrorsPropagationAnalyzer::massResolution(), and massResolution().

◆ 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

Definition at line 222 of file MuScleFitUtils.h.

Referenced by python.rootplot.utilities.Hist::__init_TGraph(), python.rootplot.utilities.Hist::__len__(), python.rootplot.root2matplotlib.Hist::_prepare_xaxis(), python.rootplot.root2matplotlib.Hist::_prepare_yaxis(), python.rootplot.utilities.Hist::av_xerr(), python.rootplot.utilities.Hist::av_yerr(), python.rootplot.utilities.Hist::delete_bin(), python.rootplot.utilities.Hist::divide(), probability(), MuScleFitBase::readProbabilityDistributionsFromFile(), python.rootplot.utilities.Hist::TGraph(), and python.rootplot.utilities.Hist::TH1F().

◆ 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

Definition at line 201 of file MuScleFitUtils.h.

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

◆ parBgrFix

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

Definition at line 205 of file MuScleFitUtils.h.

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

◆ parBgrOrder

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

Definition at line 209 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 200 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 208 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 192 of file MuScleFitUtils.h.

Referenced by ResolutionAnalyzer::analyze(), applyScale(), MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), findBestRecoRes(), massProb(), minimizeLikelihood(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().

◆ parResolFix

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

Definition at line 202 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 206 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 196 of file MuScleFitUtils.h.

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

◆ parScaleFix

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

Definition at line 203 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 207 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 190 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 251 of file MuScleFitUtils.h.

Referenced by massProb(), MuScleFit::MuScleFit(), and MuScleFitBase::readProbabilityDistributionsFromFile().

◆ 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

Definition at line 210 of file MuScleFitUtils.h.

Referenced by MuScleFit::checkParameters(), computeWeight(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), MuScleFitBase::fillHistoMap(), ResolutionAnalyzer::fillHistoMap(), findBestRecoRes(), findBestSimuRes(), findGenMuFromRes(), findSimMuFromRes(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), MuScleFitGenFilter::MuScleFitGenFilter(), MuScleFitBase::readProbabilityDistributionsFromFile(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().

◆ ResFound

bool MuScleFitUtils::ResFound = false
static

Definition at line 131 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 221 of file MuScleFitUtils.h.

Referenced by massProb(), massResolution(), probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().

◆ ResMass

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

Definition at line 136 of file MuScleFitUtils.h.

Referenced by findBestRecoRes(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), and probability().

◆ ResMaxSigma

double MuScleFitUtils::ResMaxSigma
static

Definition at line 220 of file MuScleFitUtils.h.

Referenced by massResolution(), probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().

◆ ResMinMass

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

Definition at line 137 of file MuScleFitUtils.h.

Referenced by probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().

◆ ResolFitType

int MuScleFitUtils::ResolFitType = 0
static

Definition at line 152 of file MuScleFitUtils.h.

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

◆ resolutionFunction

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

Definition at line 153 of file MuScleFitUtils.h.

Referenced by ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::massResolution(), massResolution(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().

◆ resolutionFunctionForVec

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

Definition at line 154 of file MuScleFitUtils.h.

Referenced by ResolutionAnalyzer::analyze(), MuScleFit::duringFastLoop(), minimizeLikelihood(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().

◆ 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

Definition at line 231 of file MuScleFitUtils.h.

Referenced by TestCorrection::analyze(), MuScleFit::duringFastLoop(), MuScleFit::endOfLoop(), likelihood(), minimizeLikelihood(), MuScleFit::selectMuons(), and MuScleFit::~MuScleFit().

◆ 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

Definition at line 147 of file MuScleFitUtils.h.

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

◆ speedup

bool MuScleFitUtils::speedup = false
static

Definition at line 212 of file MuScleFitUtils.h.

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

◆ 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

Definition at line 255 of file MuScleFitUtils.h.

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

◆ x

double MuScleFitUtils::x
static

Definition at line 213 of file MuScleFitUtils.h.

Referenced by python.rootplot.utilities.Hist::__init_TGraph(), svgfig.Curve.Sample::__repr__(), svgfig.Ellipse::__repr__(), python.rootplot.root2matplotlib.Hist::_prepare_xaxis(), python.rootplot.root2matplotlib.Hist::_prepare_yaxis(), applySmearing(), python.rootplot.root2matplotlib.Hist2D::box(), python.rootplot.root2matplotlib.Hist2D::contour(), python.rootplot.utilities.Hist::delete_bin(), python.rootplot.root2matplotlib.Hist::errorbar(), python.rootplot.root2matplotlib.Hist::errorbarh(), fitMass(), fitReso(), massProb(), MuScleFit::MuScleFit(), geometryXMLparser.Alignable::pos(), ntupleDataFormat._HitObject::r(), ntupleDataFormat._HitObject::r3D(), and python.rootplot.utilities.Hist::TGraph().

MuScleFitUtils::x
static double x[7][10000]
Definition: MuScleFitUtils.h:213
svgfig.canvas
def canvas(*sub, **attr)
Definition: svgfig.py:482
PDWG_BPHSkim_cff.muons
muons
Definition: PDWG_BPHSkim_cff.py:47
MuScleFitUtils::parSmear
static std::vector< double > parSmear
Definition: MuScleFitUtils.h:190
lorentzianPeak
Double_t lorentzianPeak(Double_t *x, Double_t *par)
Definition: MuScleFitUtils.cc:56
MuScleFitUtils::parScaleFix
static std::vector< int > parScaleFix
Definition: MuScleFitUtils.h:203
DDAxes::y
MuScleFitUtils::minMuonPt_
static double minMuonPt_
Definition: MuScleFitUtils.h:259
cscDigiValidation_cfi.simTrack
simTrack
Definition: cscDigiValidation_cfi.py:29
mps_fire.i
i
Definition: mps_fire.py:428
MuScleFitUtils::smearFunction
static smearFunctionBase * smearFunction
Definition: MuScleFitUtils.h:148
MuScleFitUtils::parResolMin
static std::vector< double > parResolMin
Definition: MuScleFitUtils.h:194
genParticles2HepMC_cfi.genParticles
genParticles
Definition: genParticles2HepMC_cfi.py:4
resolutionFunctionBase::covPt1Pt2
virtual double covPt1Pt2(const double &pt1, const double &eta1, const double &pt2, const double &eta2, const T &parval)
Definition: Functions.h:834
MuScleFitUtils::parvalue
static std::vector< std::vector< double > > parvalue
Definition: MuScleFitUtils.h:226
MuScleFitUtils::parScaleOrder
static std::vector< int > parScaleOrder
Definition: MuScleFitUtils.h:207
MuScleFitUtils::parorder
static std::vector< int > parorder
Definition: MuScleFitUtils.h:229
MuScleFitUtils::ResMass
static double ResMass[6]
Definition: MuScleFitUtils.h:136
muon
Definition: MuonCocktails.h:17
MuScleFitUtils::BgrFitType
static int BgrFitType
Definition: MuScleFitUtils.h:158
MuScleFitUtils::duringMinos_
static bool duringMinos_
Definition: MuScleFitUtils.h:188
MuScleFitUtils::maxMuonEtaFirstRange_
static double maxMuonEtaFirstRange_
Definition: MuScleFitUtils.h:262
resolutionFunctionBase::parNum
virtual int parNum() const
Definition: Functions.h:865
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
MuScleFitUtils::deltaPhi
static double deltaPhi(const double &phi1, const double &phi2)
Definition: MuScleFitUtils.h:109
MuScleFitUtils::massResolComponentsStruct::dmdphi1
double dmdphi1
Definition: MuScleFitUtils.h:272
MuScleFitUtils::massResolComponentsStruct::dmdcotgth2
double dmdcotgth2
Definition: MuScleFitUtils.h:275
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
mps_merge.weight
weight
Definition: mps_merge.py:88
MuScleFitUtils::parBgrFix
static std::vector< int > parBgrFix
Definition: MuScleFitUtils.h:205
gather_cfg.cout
cout
Definition: gather_cfg.py:144
MuScleFitUtils::computeMinosErrors_
static bool computeMinosErrors_
Definition: MuScleFitUtils.h:280
MuScleFitUtils::nbins
static int nbins
Definition: MuScleFitUtils.h:222
np
int np
Definition: AMPTWrapper.h:43
MuScleFitUtils::mMu2
static const double mMu2
Definition: MuScleFitUtils.h:139
BackgroundHandler::unlockParameter
bool unlockParameter(const std::vector< int > &resfind, const unsigned int ipar)
returns true if the parameter is to be unlocked
Definition: BackgroundHandler.cc:141
MuScleFitUtils::checkMassWindow
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.
Definition: MuScleFitUtils.cc:1157
MuScleFitUtils::ResHalfWidth
static double ResHalfWidth[6]
Definition: MuScleFitUtils.h:221
MuScleFitUtils::crossSectionHandler
static CrossSectionHandler * crossSectionHandler
Definition: MuScleFitUtils.h:171
MuScleFitUtils::ResMinMass
static double ResMinMass[6]
Definition: MuScleFitUtils.h:137
Gaussian
Double_t Gaussian(Double_t *x, Double_t *par)
Definition: MuScleFitUtils.cc:63
timingPdfMaker.histo
histo
Definition: timingPdfMaker.py:278
generateEDF.cont
cont
load Luminosity info ##
Definition: generateEDF.py:628
TrackCandidateProducer_cfi.simTracks
simTracks
Definition: TrackCandidateProducer_cfi.py:15
smearFunctionBase::smear
virtual void smear(double &pt, double &eta, double &phi, const double *y, const std::vector< double > &parSmear)=0
bookConverter.results
results
Definition: bookConverter.py:144
MuScleFitUtils::parResolOrder
static std::vector< int > parResolOrder
Definition: MuScleFitUtils.h:206
BackgroundHandler::setParameters
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.
Definition: BackgroundHandler.cc:102
MuScleFitUtils::parCrossSectionOrder
static std::vector< int > parCrossSectionOrder
Definition: MuScleFitUtils.h:208
resolutionFunctionBase::sigmaCotgTh
virtual double sigmaCotgTh(const double &pt, const double &eta, const T &parval)=0
MuScleFitUtils::applyScale
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)
Definition: MuScleFitUtils.cc:465
MuScleFitUtils::likelihoodInLoop_
static TH1D * likelihoodInLoop_
Definition: MuScleFitUtils.h:184
ResolutionFunction::sigmaPhi
double sigmaPhi(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
Definition: ResolutionFunction.h:85
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
MuScleFitUtils::FitStrategy
static int FitStrategy
Definition: MuScleFitUtils.h:211
MuScleFitUtils::massResolComponents
static struct MuScleFitUtils::massResolComponentsStruct massResolComponents
Definition: MuScleFitUtils.cc:262
likelihood
void likelihood(int &npar, double *grad, double &fval, double *xval, int flag)
Definition: MuScleFitUtils.cc:1784
MuScleFitUtils::useProbsFile_
static bool useProbsFile_
Definition: MuScleFitUtils.h:255
MuScleFitUtils::resfind
static std::vector< int > resfind
Definition: MuScleFitUtils.h:210
GL
TF1 * GL
Definition: MuScleFitUtils.cc:83
chg
const float chg[109]
Definition: CoreSimTrack.cc:5
HepMC::GenEvent
Definition: hepmc_rootio.cc:9
MuScleFitUtils::rapidityBinsForZ_
static bool rapidityBinsForZ_
Definition: MuScleFitUtils.h:251
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MuScleFitUtils::normalizationChanged_
static unsigned int normalizationChanged_
Definition: MuScleFitUtils.h:245
MuScleFitUtils::minimumShapePlots_
static bool minimumShapePlots_
Definition: MuScleFitUtils.h:281
MuScleFitUtils::doCrossSectionFit
static std::vector< int > doCrossSectionFit
Definition: MuScleFitUtils.h:180
ires
int ires[2]
Definition: CascadeWrapper.h:19
MuScleFitUtils::massResolComponentsStruct::dmdpt2
double dmdpt2
Definition: MuScleFitUtils.h:271
RPCNoise_example.check
check
Definition: RPCNoise_example.py:71
part
part
Definition: HCALResponse.h:20
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
BackgroundHandler::regionsParNum
int regionsParNum()
Returns the total number of parameters used for the regions.
Definition: BackgroundHandler.h:48
MuScleFitUtils::parResol
static std::vector< double > parResol
Definition: MuScleFitUtils.h:192
MuScleFitUtils::parCrossSection
static std::vector< double > parCrossSection
Definition: MuScleFitUtils.h:200
MuScleFitUtils::debugMassResol_
static bool debugMassResol_
Definition: MuScleFitUtils.h:268
PVValHelper::eta
Definition: PVValidationHelpers.h:70
MuScleFitUtils::parResolMax
static std::vector< double > parResolMax
Definition: MuScleFitUtils.h:195
w
const double w
Definition: UKUtility.cc:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HLT_FULL_cff.eta2
eta2
Definition: HLT_FULL_cff.py:9551
CrossSectionHandler::setParameters
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.
Definition: CrossSectionHandler.h:58
MuScleFitUtils::ResFound
static bool ResFound
Definition: MuScleFitUtils.h:131
MuScleFitUtils::motherPdgIdArray
static const unsigned int motherPdgIdArray[6]
Definition: MuScleFitUtils.h:143
HLT_FULL_cff.pt1
pt1
Definition: HLT_FULL_cff.py:9893
MuScleFitUtils::invDimuonMass
static double invDimuonMass(const lorentzVector &mu1, const lorentzVector &mu2)
Definition: MuScleFitUtils.cc:518
MuScleFitUtils::scaleFunction
static scaleFunctionBase< double * > * scaleFunction
Definition: MuScleFitUtils.h:156
MuScleFitUtils::massResolution
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
MuScleFitUtils::startWithSimplex_
static bool startWithSimplex_
Definition: MuScleFitUtils.h:279
MuScleFitUtils::minMuonEtaSecondRange_
static double minMuonEtaSecondRange_
Definition: MuScleFitUtils.h:263
MuScleFitUtils::ResolFitType
static int ResolFitType
Definition: MuScleFitUtils.h:152
MuScleFitUtils::parScale
static std::vector< double > parScale
Definition: MuScleFitUtils.h:196
MuScleFitUtils::resolutionFunctionForVec
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
Definition: MuScleFitUtils.h:154
MuScleFitUtils::minMuonEtaFirstRange_
static double minMuonEtaFirstRange_
Definition: MuScleFitUtils.h:261
HLT_FULL_cff.eta1
eta1
Definition: HLT_FULL_cff.py:9550
MuScleFitUtils::massWindowHalfWidth
static double massWindowHalfWidth[3][6]
Definition: MuScleFitUtils.h:134
resolutionFunctionBase::sigmaPhi
virtual double sigmaPhi(const double &pt, const double &eta, const T &parval)=0
runTauDisplay.gp
gp
Definition: runTauDisplay.py:431
MuScleFitUtils::massResolComponentsStruct::dmdcotgth1
double dmdcotgth1
Definition: MuScleFitUtils.h:274
scaleFunctionBase::setParameters
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.
MuScleFitUtils::debug
static int debug
Definition: MuScleFitUtils.h:130
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
MuScleFitUtils::parScaleMax
static std::vector< double > parScaleMax
Definition: MuScleFitUtils.h:199
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
lorentzVector
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:9
BackgroundHandler::backgroundFunction
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)
Definition: BackgroundHandler.cc:249
iorder
int iorder
Definition: CascadeWrapper.h:49
scaleFunctionBase::parNum
virtual int parNum() const
Definition: Functions.h:76
MuScleFitUtils::parfix
static std::vector< int > parfix
Definition: MuScleFitUtils.h:228
MuScleFitUtils::SavedPair
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
Definition: MuScleFitUtils.h:231
createfilelist.int
int
Definition: createfilelist.py:10
scaleFunctionBase::resetParameters
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
resolutionFunctionBase::sigmaPt
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
ResolutionFunction::sigmaPt
double sigmaPt(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
Definition: ResolutionFunction.h:65
MuScleFitUtils::ResGamma
static double ResGamma[6]
Definition: MuScleFitUtils.h:135
MuScleFitUtils::biasFunction
static scaleFunctionBase< std::vector< double > > * biasFunction
Definition: MuScleFitUtils.h:151
scaleFunctionBase::scale
virtual double scale(const double &pt, const double &eta, const double &phi, const int chg, const T &parScale) const =0
HLT_FULL_cff.pt2
pt2
Definition: HLT_FULL_cff.py:9895
edm::HepMCProduct::GetEvent
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
MuScleFitUtils::GLNorm
static double GLNorm[6][1001]
Definition: MuScleFitUtils.h:219
MuScleFitUtils::parResolFix
static std::vector< int > parResolFix
Definition: MuScleFitUtils.h:202
submitPVResolutionJobs.err
err
Definition: submitPVResolutionJobs.py:85
resolutionFunctionBase::setParameters
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
MuScleFitUtils::GLZNorm
static double GLZNorm[40][1001]
Definition: MuScleFitUtils.h:217
MuScleFitUtils::parBgr
static std::vector< double > parBgr
Definition: MuScleFitUtils.h:201
MuScleFitUtils::probability
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...
Definition: MuScleFitUtils.cc:750
MuScleFitUtils::maxMuonPt_
static double maxMuonPt_
Definition: MuScleFitUtils.h:260
MuScleFitUtils::rminPtr_
static TMinuit * rminPtr_
Definition: MuScleFitUtils.h:242
MuScleFitUtils::parBgrOrder
static std::vector< int > parBgrOrder
Definition: MuScleFitUtils.h:209
CrossSectionHandler::relativeCrossSections
std::vector< double > relativeCrossSections(const double *variables, const std::vector< int > &resfind)
Perform a variable transformation from N-1 to relative cross sections.
Definition: CrossSectionHandler.h:124
groupFilesInBlocks.nn
nn
Definition: groupFilesInBlocks.py:150
DDAxes::phi
MuScleFitUtils::parBias
static std::vector< double > parBias
Definition: MuScleFitUtils.h:191
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
BackgroundHandler::windowBorders
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...
Definition: BackgroundHandler.cc:170
MuScleFitUtils::scaleFunctionForVec
static scaleFunctionBase< std::vector< double > > * scaleFunctionForVec
Definition: MuScleFitUtils.h:157
MuScleFitUtils::ResMaxSigma
static double ResMaxSigma[6]
Definition: MuScleFitUtils.h:220
GenParticle.GenParticle
GenParticle
Definition: GenParticle.py:18
MuScleFitUtils::doBackgroundFit
static std::vector< int > doBackgroundFit
Definition: MuScleFitUtils.h:181
CrossSectionHandler::parNum
unsigned int parNum()
Definition: CrossSectionHandler.h:121
MuScleFitUtils::minuitLoop_
static int minuitLoop_
Definition: MuScleFitUtils.h:183
edm::shift
static unsigned const int shift
Definition: LuminosityBlockID.cc:7
MuScleFitUtils::ScaleFitType
static int ScaleFitType
Definition: MuScleFitUtils.h:155
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:29
MuScleFitUtils::signalProb_
static TH1D * signalProb_
Definition: MuScleFitUtils.h:185
MuScleFitUtils::muMass
static const double muMass
Definition: MuScleFitUtils.h:140
MuScleFitUtils::GLZValue
static double GLZValue[40][1001][1001]
Definition: MuScleFitUtils.h:216
MuScleFitUtils::backgroundProb_
static TH1D * backgroundProb_
Definition: MuScleFitUtils.h:186
MuScleFitUtils::SmearType
static int SmearType
Definition: MuScleFitUtils.h:147
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
MuScleFitUtils::GLValue
static double GLValue[6][1001][1001]
Definition: MuScleFitUtils.h:218
MuScleFitUtils::counter_resprob
static int counter_resprob
Definition: MuScleFitUtils.h:215
MuScleFitUtils::fromPtEtaPhiToPxPyPz
static lorentzVector fromPtEtaPhiToPxPyPz(const double *ptEtaPhiE)
Definition: MuScleFitUtils.cc:502
BackgroundHandler::rescale
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.)
Definition: BackgroundHandler.cc:192
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
MuScleFitUtils::parScaleMin
static std::vector< double > parScaleMin
Definition: MuScleFitUtils.h:198
MuScleFitUtils::massResolComponentsStruct::dmdphi2
double dmdphi2
Definition: MuScleFitUtils.h:273
MuScleFitUtils::scaleFitNotDone_
static bool scaleFitNotDone_
Definition: MuScleFitUtils.h:238
BeamSpotPI::Y
Definition: BeamSpotPayloadInspectorHelper.h:32
MuScleFitUtils::doScaleFit
static std::vector< int > doScaleFit
Definition: MuScleFitUtils.h:179
BackgroundHandler::resMass
double resMass(const bool doBackgroundFit, const int ires)
Definition: BackgroundHandler.cc:182
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
MuScleFitUtils::parScaleStep
static std::vector< double > parScaleStep
Definition: MuScleFitUtils.h:197
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
MuScleFitUtils::loopCounter
static unsigned int loopCounter
Definition: MuScleFitUtils.h:145
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
P
std::pair< OmniClusterRef, TrackingParticleRef > P
Definition: BDHadronTrackMonitoringAnalyzer.cc:203
trackingParticleMuon_cfi.pmin
pmin
Definition: trackingParticleMuon_cfi.py:10
CrossSectionHandler::releaseParameters
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.
Definition: CrossSectionHandler.h:97
MuScleFitUtils::sherpa_
static bool sherpa_
Definition: MuScleFitUtils.h:248
ResolutionFunction::sigmaCotgTh
double sigmaCotgTh(const U &track, const int i=0) const
The second, optional, parameter is the iteration number.
Definition: ResolutionFunction.h:75
MuScleFitUtils::parResolStep
static std::vector< double > parResolStep
Definition: MuScleFitUtils.h:193
MuScleFitUtils::resolutionFunction
static resolutionFunctionBase< double * > * resolutionFunction
Definition: MuScleFitUtils.h:153
MuScleFitUtils::massResolComponentsStruct::dmdpt1
double dmdpt1
Definition: MuScleFitUtils.h:270
MuScleFitUtils::MuonType
static int MuonType
Definition: MuScleFitUtils.h:223
TtFullHadEvtBuilder_cfi.prob
prob
Definition: TtFullHadEvtBuilder_cfi.py:33
geometryCSVtoXML.xx
xx
Definition: geometryCSVtoXML.py:19
weight
Definition: weight.py:1
MuScleFitUtils::goodmuon
static int goodmuon
Definition: MuScleFitUtils.h:214
MuScleFitUtils::backgroundHandler
static BackgroundHandler * backgroundHandler
Definition: MuScleFitUtils.h:175
MuScleFitUtils::doResolFit
static std::vector< int > doResolFit
Definition: MuScleFitUtils.h:178
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
MuScleFitUtils::ReducedSavedPair
static std::vector< std::pair< lorentzVector, lorentzVector > > ReducedSavedPair
Definition: MuScleFitUtils.h:232
MuScleFitUtils::massProb
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
Definition: MuScleFitUtils.cc:2208
MuScleFitUtils::totalResNum
static const int totalResNum
Definition: MuScleFitUtils.h:133
Generated for CMSSW Reference Manual by doxygen 1.8.16