#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.
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.
static double	deltaR (const double &eta1, const double &eta2, const double &phi1, const double &phi2)
static std::pair < lorentzVector, lorentzVector >	findBestRecoRes (const std::vector< reco::LeafCandidate > &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=false)
static double	massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, double *parval, const bool doUseBkgrWindow=false)
static double	massResolution (const lorentzVector &mu1, const lorentzVector &mu2, double *parval)
static double	massResolution (const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval)
static double	massResolution (const lorentzVector &mu1, const lorentzVector &mu2)
static double	massResolution (const lorentzVector &mu1, const lorentzVector &mu2, std::auto_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.
Static Public Attributes
static const int	backgroundFunctionsRegions
static BackgroundHandler *	backgroundHandler
static TH1D *	backgroundProb_ = 0
static int	BgrFitType = 0
static scaleFunctionBase < std::vector< double > > *	biasFunction = 0
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 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 < 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_ = 0
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< int >	parResolOrder
static std::vector< double >	parScale
static std::vector< int >	parScaleFix
static std::vector< int >	parScaleOrder
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 = 0
static resolutionFunctionBase < std::vector< double > > *	resolutionFunctionForVec = 0
static TMinuit *	rminPtr_ = 0
static std::vector< std::pair < lorentzVector, lorentzVector > >	SavedPair
static bool	scaleFitNotDone_ = true
static int	ScaleFitType = 0
static scaleFunctionBase < double * > *	scaleFunction = 0
static scaleFunctionBase < std::vector< double > > *	scaleFunctionForVec = 0
static bool	sherpa_ = false
static TH1D *	signalProb_ = 0
static std::vector< std::pair < lorentzVector, lorentzVector > >	simPair
static smearFunctionBase *	smearFunction = 0
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 46 of file MuScleFitUtils.h.

Constructor & Destructor Documentation

MuScleFitUtils::MuScleFitUtils ( ) [inline]

Definition at line 51 of file MuScleFitUtils.h.

{};

virtual MuScleFitUtils::~MuScleFitUtils ( ) [inline, virtual]

Definition at line 55 of file MuScleFitUtils.h.

{};

Member Function Documentation

lorentzVector MuScleFitUtils::applyBias	(	const lorentzVector &	muon,
		const int	charge
	)		`[static]`

Definition at line 417 of file MuScleFitUtils.cc.

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

{
  double ptEtaPhiE[4] = {muon.Pt(),muon.Eta(),muon.Phi(),muon.E()};

  if (MuScleFitUtils::debug>1) std::cout << "pt before bias = " << ptEtaPhiE[0] << std::endl;

  // Use functors (although not with the () operator)
  // Note that we always pass pt, eta and phi, but internally only the needed
  // values are used.
  // The functors used are takend from the same group used for the scaling
  // thus the name of the method used is "scale".
  ptEtaPhiE[0] = biasFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, MuScleFitUtils::parBias);

  if (MuScleFitUtils::debug>1) std::cout << "pt after bias = " << ptEtaPhiE[0] << std::endl;

  return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}

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

Definition at line 437 of file MuScleFitUtils.cc.

References ExpressReco_HICollisions_FallBack::id, and L1TEmulatorMonitor_cff::p.

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

{
  double * p = new double[(int)(parval.size())];
  // Replaced by auto_ptr, which handles delete at the end
  // std::auto_ptr<double> p(new double[(int)(parval.size())]);
  // Removed auto_ptr, check massResolution for an explanation.
  int id = 0;
  for (std::vector<double>::const_iterator it=parval.begin(); it!=parval.end(); ++it, ++id) {
    //(&*p)[id] = *it;
    // Also ok would be (p.get())[id] = *it;
    p[id] = *it;
  }
  lorentzVector tempScaleVec( applyScale (muon, p, chg) );
  delete[] p;
  return tempScaleVec;
}

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

Definition at line 457 of file MuScleFitUtils.cc.

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

{
  double ptEtaPhiE[4] = {muon.Pt(),muon.Eta(),muon.Phi(),muon.E()};
  int shift = parResol.size();

  if (MuScleFitUtils::debug>1) std::cout << "pt before scale = " << ptEtaPhiE[0] << std::endl;

  // the address of parval[shift] is passed as pointer to double. Internally it is used as a normal array, thus:
  // array[0] = parval[shift], array[1] = parval[shift+1], ...
  ptEtaPhiE[0] = scaleFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, &(parval[shift]));

  if (MuScleFitUtils::debug>1) std::cout << "pt after scale = " << ptEtaPhiE[0] << std::endl;

  return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}

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

Definition at line 387 of file MuScleFitUtils.cc.

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

{
  double pt = muon.Pt();
  double eta = muon.Eta();
  double phi = muon.Phi();
  double E = muon.E();

  double y[7];
  for (int i=0; i<SmearType+3; i++) {
    y[i] = x[i][goodmuon%10000];
  }

  // Use the smear function selected in the constructor
  smearFunction->smear( pt, eta, phi, y, parSmear );

  if (debug>9) {
    std::cout << "Smearing Pt,eta,phi = " << pt << " " <<  eta << " "
         << phi << "; x = ";
    for (int i=0; i<SmearType+3; i++) {
      std::cout << y[i];
    }
    std::cout << std::endl;
  }

  double ptEtaPhiE[4] = {pt, eta, phi, E};
  return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}

bool MuScleFitUtils::checkMassWindow	(	const double &	mass,
		const double &	leftBorder,
		const double &	rightBorder
	)		`[inline, static]`

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

Definition at line 1072 of file MuScleFitUtils.cc.

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

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

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

Definition at line 1079 of file MuScleFitUtils.cc.

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

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

{
  // Compute weight for this event
  // -----------------------------
  double weight = 0.;

  // Take the highest-mass resonance within bounds
  // NB this must be revised once credible estimates of the relative xs of Y(1S), (2S), and (3S)
  // are made. Those are priors in the decision of which resonance to assign to an in-between event.
  // -----------------------------------------------------------------------------------------------

  if( doUseBkgrWindow && (debug > 0) ) std::cout << "using backgrond window for mass = " << mass << std::endl;
  // bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
  bool useBackgroundWindow = (doBackgroundFit[loopCounter]);

  for (int ires=0; ires<6; ires++) {
    if (resfind[ires]>0 && weight==0.) {
      // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
      std::pair<double, double> windowBorder = backgroundHandler->windowBorders( useBackgroundWindow, ires );
      // if( checkMassWindow(mass, ires, backgroundHandler->resMass( useBackgroundWindow, ires ),
      //                     windowFactor.first, windowFactor.second) ) {
      if( checkMassWindow(mass, windowBorder.first, windowBorder.second) ) {
        weight = 1.0;
        if( doUseBkgrWindow && (debug > 0) ) std::cout << "setting weight to = " << weight << std::endl;
      }
    }
  }

  return weight;
}

static double MuScleFitUtils::deltaPhi	(	const double &	phi1,
		const double &	phi2
	)		`[inline, static]`

Definition at line 88 of file MuScleFitUtils.h.

References Pi.

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

  {
    double deltaPhi = phi1 - phi2;
    while(deltaPhi >= TMath::Pi()) deltaPhi -= 2*TMath::Pi();
    while(deltaPhi < -TMath::Pi()) deltaPhi += 2*TMath::Pi();
    return fabs(deltaPhi);
  }

static double MuScleFitUtils::deltaPhiNoFabs	(	const double &	phi1,
		const double &	phi2
	)		`[inline, static]`

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

Definition at line 96 of file MuScleFitUtils.h.

References deltaPhi(), and Pi.

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

  {
    double deltaPhi = phi1 - phi2;
    while(deltaPhi >= TMath::Pi()) deltaPhi -= 2*TMath::Pi();
    while(deltaPhi < -TMath::Pi()) deltaPhi += 2*TMath::Pi();
    return deltaPhi;
  }

static double MuScleFitUtils::deltaR	(	const double &	eta1,
		const double &	eta2,
		const double &	phi1,
		const double &	phi2
	)		`[inline, static]`

Definition at line 103 of file MuScleFitUtils.h.

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

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

std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findBestRecoRes ( const std::vector< reco::LeafCandidate > & muons ) [static]

Definition at line 297 of file MuScleFitUtils.cc.

References gather_cfg::cout, debug, massProb(), massResolution(), maxMuonEtaFirstRange_, maxMuonEtaSecondRange_, maxMuonPt_, minMuonEtaFirstRange_, minMuonEtaSecondRange_, minMuonPt_, parResol, resfind, ResFound, ResMass, and useProbsFile_.

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

                                                                                                               {
  // NB this routine returns the resonance, but it also sets the ResFound flag, which
  // is used in MuScleFit to decide whether to use the event or not.
  // --------------------------------------------------------------------------------
  if (debug>0) std::cout << "In findBestRecoRes" << std::endl;
  ResFound = false;
  std::pair<lorentzVector, lorentzVector> recMuFromBestRes;

  // Choose the best resonance using its mass probability
  // ----------------------------------------------------
  double maxprob = -0.1;
  double minDeltaMass = 999999;
  std::pair<reco::LeafCandidate,reco::LeafCandidate> bestMassMuons;
  for (std::vector<reco::LeafCandidate>::const_iterator Muon1=muons.begin(); Muon1!=muons.end(); ++Muon1) {
    for (std::vector<reco::LeafCandidate>::const_iterator Muon2=Muon1+1; Muon2!=muons.end(); ++Muon2) {
      if (((*Muon1).charge()*(*Muon2).charge())>0) {
        continue; // This also gets rid of Muon1==Muon2...
      }
      // To allow the selection of ranges at negative and positive eta independently we define two
      // ranges of eta: (minMuonEtaFirstRange_, maxMuonEtaFirstRange_) and (minMuonEtaSecondRange_, maxMuonEtaSecondRange_).
      // If the interval selected is simmetric, one only needs to specify the first range. The second has
      // default values that accept all muons (minMuonEtaSecondRange_ = -100., maxMuonEtaSecondRange_ = 100.).
      double pt1 = (*Muon1).p4().Pt();
      double pt2 = (*Muon2).p4().Pt();
      double eta1 = (*Muon1).p4().Eta();
      double eta2 = (*Muon2).p4().Eta();
      if( pt1 > minMuonPt_ && pt1 < maxMuonPt_ &&
          pt2 > minMuonPt_ && pt2 < maxMuonPt_ &&
          eta1 > minMuonEtaFirstRange_ && eta1 < maxMuonEtaFirstRange_ &&
          eta2 > minMuonEtaFirstRange_ && eta2 < maxMuonEtaFirstRange_ &&
          eta1 > minMuonEtaSecondRange_ && eta1 < maxMuonEtaSecondRange_ &&
          eta2 > minMuonEtaSecondRange_ && eta2 < maxMuonEtaSecondRange_ ) {
        double mcomb = ((*Muon1).p4()+(*Muon2).p4()).mass();
        double Y = ((*Muon1).p4()+(*Muon2).p4()).Rapidity();
        if (debug>1) {
          std::cout<<"muon1 "<<(*Muon1).p4().Px()<<", "<<(*Muon1).p4().Py()<<", "<<(*Muon1).p4().Pz()<<", "<<(*Muon1).p4().E()<<std::endl;
          std::cout<<"muon2 "<<(*Muon2).p4().Px()<<", "<<(*Muon2).p4().Py()<<", "<<(*Muon2).p4().Pz()<<", "<<(*Muon2).p4().E()<<std::endl;
          std::cout<<"mcomb "<<mcomb<<std::endl;}
        double massResol = 0.;
        if( useProbsFile_ ) {
          massResol = massResolution ((*Muon1).p4(), (*Muon2).p4(), parResol);
        }
        double prob = 0;
        for( int ires=0; ires<6; ires++ ) {
          if( resfind[ires]>0 ) {
            if( useProbsFile_ ) {
              prob = massProb( mcomb, Y, ires, massResol );
            }
            if( prob>maxprob ) {
              if( (*Muon1).charge()<0 ) { // store first the mu minus and then the mu plus
                recMuFromBestRes.first = (*Muon1).p4();
                recMuFromBestRes.second = (*Muon2).p4();
              } else {
                recMuFromBestRes.first = (*Muon2).p4();
                recMuFromBestRes.second = (*Muon1).p4();
              }
              ResFound = true; // NNBB we accept "resonances" even outside mass bounds
              maxprob = prob;
            }
            // if( ResMass[ires] == 0 ) {
            //   std::cout << "Error: ResMass["<<ires<<"] = " << ResMass[ires] << std::endl;
            //   exit(1);
            // }
            double deltaMass = fabs(mcomb-ResMass[ires])/ResMass[ires];
            if( deltaMass<minDeltaMass ){
              bestMassMuons = std::make_pair((*Muon1),(*Muon2));
              minDeltaMass = deltaMass;
            }
          }
        }
      }
    }
  }
  //If outside mass window (maxprob==0) then take the two muons with best invariant mass
  //(anyway they will not be used in the likelihood calculation, only to fill plots)
  if(!maxprob){
    if(bestMassMuons.first.charge()<0){
      recMuFromBestRes.first = bestMassMuons.first.p4();
      recMuFromBestRes.second = bestMassMuons.second.p4();
    }
    else{
      recMuFromBestRes.second = bestMassMuons.first.p4();
      recMuFromBestRes.first = bestMassMuons.second.p4();
    }
  }
  return recMuFromBestRes;
}

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

Definition at line 260 of file MuScleFitUtils.cc.

References massProb(), and resfind.

Referenced by MuScleFitPlotter::fillSim().

                                                                                               {

  std::pair<SimTrack, SimTrack> simMuFromBestRes;
  double maxprob = -0.1;

  // Double loop on muons
  // --------------------
  for (std::vector<SimTrack>::const_iterator simMu1=simMuons.begin(); simMu1!=simMuons.end(); simMu1++) {
    for (std::vector<SimTrack>::const_iterator simMu2=simMu1+1; simMu2!=simMuons.end(); simMu2++) {
      if (((*simMu1).charge()*(*simMu2).charge())>0) {
        continue; // this also gets rid of simMu1==simMu2...
      }
      // Choose the best resonance using its mass. Check Z, Y(3S,2S,1S), Psi(2S), J/Psi in order
      // ---------------------------------------------------------------------------------------
      double mcomb = ((*simMu1).momentum()+(*simMu2).momentum()).mass();
      double Y = ((*simMu1).momentum()+(*simMu2).momentum()).Rapidity();
      for (int ires=0; ires<6; ires++) {
        if (resfind[ires]>0) {
          double prob = massProb( mcomb, Y, ires, 0. );
          if (prob>maxprob) {
            simMuFromBestRes.first = (*simMu1);
            simMuFromBestRes.second = (*simMu2);
            maxprob = prob;
          }
        }
      }
    }
  }

  // Return most likely combination of muons making a resonance
  // ----------------------------------------------------------
  return simMuFromBestRes;
}

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

Definition at line 2164 of file MuScleFitUtils.cc.

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

{
  const HepMC::GenEvent* Evt = evtMC->GetEvent();
  std::pair<lorentzVector,lorentzVector> muFromRes;
  //Loop on generated particles
  for (HepMC::GenEvent::particle_const_iterator part=Evt->particles_begin();
       part!=Evt->particles_end(); part++) {
    if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
      bool fromRes = false;
      for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);
           mother != (*part)->production_vertex()->particles_end(HepMC::ancestors); ++mother) {
        unsigned int motherPdgId = (*mother)->pdg_id();

        // For sherpa the resonance is not saved. The muons from the resonance can be identified
        // by having as mother a muon of status 3.
        if( sherpa_ ) {
          if( motherPdgId == 13 && (*mother)->status() == 3 ) fromRes = true;
        }
        else {
          for( int ires = 0; ires < 6; ++ires ) {
            if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
          }
        }
      }
      if(fromRes){
        if((*part)->pdg_id()==13)
          //   muFromRes.first = (*part)->momentum();
          muFromRes.first = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
                                           (*part)->momentum().pz(),(*part)->momentum().e()));
        else
          muFromRes.second = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
                                            (*part)->momentum().pz(),(*part)->momentum().e()));
      }
    }
  }
  return muFromRes;
}

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

Definition at line 2202 of file MuScleFitUtils.cc.

References gather_cfg::cout, debug, motherPdgIdArray, and resfind.

Referenced by MuScleFitGenFilter::filter().

{
  std::pair<lorentzVector,lorentzVector> muFromRes;

  //Loop on generated particles
  if( debug>0 ) std::cout << "Starting loop on " << genParticles->size() << " genParticles" << std::endl;
  for( reco::GenParticleCollection::const_iterator part=genParticles->begin(); part!=genParticles->end(); ++part ) {
    if (fabs(part->pdgId())==13 && part->status()==1) {
      bool fromRes = false;
      unsigned int motherPdgId = part->mother()->pdgId();
      if( debug>0 ) {
        std::cout << "Found a muon with mother: " << motherPdgId << std::endl;
      }
      for( int ires = 0; ires < 6; ++ires ) {
        if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
      }
      if(fromRes){
        if(part->pdgId()==13) {
          muFromRes.first = part->p4();
          if( debug>0 ) std::cout << "Found a genMuon + : " << muFromRes.first << std::endl;
          //      muFromRes.first = (lorentzVector(part->p4().px(),part->p4().py(),
          //                                       part->p4().pz(),part->p4().e()));
        }
        else {
          muFromRes.second = part->p4();
          if( debug>0 ) std::cout << "Found a genMuon - : " << muFromRes.second << std::endl;
          //      muFromRes.second = (lorentzVector(part->p4().px(),part->p4().py(),
          //                                        part->p4().pz(),part->p4().e()));
        }
      }
    }
  }
  return muFromRes;
}

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

Definition at line 2126 of file MuScleFitUtils.cc.

References configurableAnalysis::GenParticle, motherPdgIdArray, and resfind.

Referenced by MuScleFitPlotter::fillGenSim().

{
  //Loop on simulated tracks
  std::pair<lorentzVector, lorentzVector> simMuFromRes;
  for( edm::SimTrackContainer::const_iterator simTrack=simTracks->begin(); simTrack!=simTracks->end(); ++simTrack ) {
    //Chose muons
    if (fabs((*simTrack).type())==13) {
      //If tracks from IP than find mother
      if ((*simTrack).genpartIndex()>0) {
        HepMC::GenParticle* gp = evtMC->GetEvent()->barcode_to_particle ((*simTrack).genpartIndex());
        if( gp != 0 ) {

          for (HepMC::GenVertex::particle_iterator mother = gp->production_vertex()->particles_begin(HepMC::ancestors);
               mother!=gp->production_vertex()->particles_end(HepMC::ancestors); ++mother) {

            bool fromRes = false;
            unsigned int motherPdgId = (*mother)->pdg_id();
            for( int ires = 0; ires < 6; ++ires ) {
              if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
            }
            if( fromRes ) {
              if(gp->pdg_id() == 13)
                simMuFromRes.first = lorentzVector(simTrack->momentum().px(),simTrack->momentum().py(),
                                                   simTrack->momentum().pz(),simTrack->momentum().e());
              else
                simMuFromRes.second = lorentzVector(simTrack->momentum().px(),simTrack->momentum().py(),
                                                    simTrack->momentum().pz(),simTrack->momentum().e());
            }
          }
        }
        // else LogDebug("MuScleFitUtils") << "WARNING: no matching genParticle found for simTrack" << std::endl;
      }
    }
  }
  return simMuFromRes;
}

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

Definition at line 1845 of file MuScleFitUtils.cc.

References ExpressReco_HICollisions_FallBack::chi2, cont, gather_cfg::cout, debug, ExpressReco_HICollisions_FallBack::e, i, j, lorentzianPeak(), AlCaRecoCosmics_cfg::name, Gflash::par, python::entryComment::results, and x.

                                                             {

  if (MuScleFitUtils::debug>0) std::cout << "Fitting " << histo->GetName() << std::endl;

  std::vector<TGraphErrors *> results;

  // Results of the fit
  // ------------------
  std::vector<double> Ftop;
  std::vector<double> Fwidth;
  std::vector<double> Fmass;
  std::vector<double> Etop;
  std::vector<double> Ewidth;
  std::vector<double> Emass;
  std::vector<double> Fchi2;
  // X bin center and width
  // ----------------------
  std::vector<double> Xcenter;
  std::vector<double> Ex;

  // Fit with lorentzian peak
  // ------------------------
  TF1 *fitFcn = new TF1 ("fitFcn", lorentzianPeak, 70, 110, 3);
  fitFcn->SetParameters (100, 3, 91);
  fitFcn->SetParNames ("Ftop", "Fwidth", "Fmass");
  fitFcn->SetLineWidth (2);

  // Fit slices projected along Y from bins in X
  // -------------------------------------------
  double cont_min = 20;    // Minimum number of entries
  Int_t binx = histo->GetXaxis()->GetNbins();
  // TFile *f= new TFile("prova.root", "recreate");
  // histo->Write();
  for (int i=1; i<=binx; i++) {
    TH1 * histoY = histo->ProjectionY ("", i, i);
    // histoY->Write();
    double cont = histoY->GetEntries();
    if (cont>cont_min) {
      histoY->Fit ("fitFcn", "0", "", 70, 110);
      double *par = fitFcn->GetParameters();
      double *err = fitFcn->GetParErrors();

      Ftop.push_back(par[0]);
      Fwidth.push_back(par[1]);
      Fmass.push_back(par[2]);
      Etop.push_back(err[0]);
      Ewidth.push_back(err[1]);
      Emass.push_back(err[2]);

      double chi2 = fitFcn->GetChisquare();
      Fchi2.push_back(chi2);

      double xx = histo->GetXaxis()->GetBinCenter(i);
      Xcenter.push_back(xx);
      double ex = 0; // FIXME: you can use the bin width
      Ex.push_back(ex);
    }
  }
  // f->Close();

  // Put the fit results in arrays for TGraphErrors
  // ----------------------------------------------
  const int nn = Fmass.size();
  double *x    = new double[nn];
  double *ym   = new double[nn];
  double *e    = new double[nn];
  double *eym  = new double[nn];
  double *yw   = new double[nn];
  double *eyw  = new double[nn];
  double *yc   = new double[nn];

  for (int j=0; j<nn; j++) {
    x[j]   = Xcenter[j];
    ym[j]  = Fmass[j];
    eym[j] = Emass[j];
    yw[j]  = Fwidth[j];
    eyw[j] = Ewidth[j];
    yc[j]  = Fchi2[j];
    e[j]   = Ex[j];
  }

  // Create TGraphErrors
  // -------------------
  TString name = histo->GetName();
  TGraphErrors *grM = new TGraphErrors (nn, x, ym, e, eym);
  grM->SetTitle (name+"_M");
  grM->SetName (name+"_M");
  TGraphErrors *grW = new TGraphErrors (nn, x, yw, e, eyw);
  grW->SetTitle (name+"_W");
  grW->SetName (name+"_W");
  TGraphErrors *grC = new TGraphErrors (nn, x, yc, e, e);
  grC->SetTitle (name+"_chi2");
  grC->SetName (name+"_chi2");

  // Cleanup
  // -------
  delete x;
  delete ym;
  delete eym;
  delete yw;
  delete eyw;
  delete yc;
  delete e;
  delete fitFcn;

  results.push_back(grM);
  results.push_back(grW);
  results.push_back(grC);
  return results;
}

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

Definition at line 1958 of file MuScleFitUtils.cc.

References ExpressReco_HICollisions_FallBack::chi2, cont, gather_cfg::cout, ExpressReco_HICollisions_FallBack::e, Gaussian(), i, j, AlCaRecoCosmics_cfg::name, Gflash::par, python::entryComment::results, and x.

                                                             {
  std::cout << "Fitting " << histo->GetName() << std::endl;
  std::vector<TGraphErrors *> results;

  // Results from fit
  // ----------------
  std::vector<double> maxs;
  std::vector<double> means;
  std::vector<double> sigmas;
  std::vector<double> chi2s;
  std::vector<double> Emaxs;
  std::vector<double> Emeans;
  std::vector<double> Esigmas;

  // X bin center and width
  // ----------------------
  std::vector<double> Xcenter;
  std::vector<double> Ex;

  // Fit with a gaussian
  // -------------------
  TF1 *fitFcn = new TF1 ("fitFunc", Gaussian, -0.2, 0.2, 3);
  fitFcn->SetParameters (100, 0, 0.02);
  fitFcn->SetParNames ("max", "mean", "sigma");
  fitFcn->SetLineWidth (2);

  // Fit slices projected along Y from bins in X
  // -------------------------------------------
  double cont_min = 20;    // Minimum number of entries
  Int_t binx = histo->GetXaxis()->GetNbins();
  for (int i=1; i<=binx; i++) {
    TH1 * histoY = histo->ProjectionY ("", i, i);
    double cont = histoY->GetEntries();
    if (cont>cont_min) {
      histoY->Fit ("fitFunc", "0", "", -0.2, 0.2);
      double *par = fitFcn->GetParameters();
      double *err = fitFcn->GetParErrors();

      maxs.push_back (par[0]);
      means.push_back (par[1]);
      sigmas.push_back (par[2]);
      Emaxs.push_back (err[0]);
      Emeans.push_back (err[1]);
      Esigmas.push_back (err[2]);

      double chi2 = fitFcn->GetChisquare();
      chi2s.push_back (chi2);

      double xx = histo->GetXaxis()->GetBinCenter(i);
      Xcenter.push_back (xx);
      double ex = 0; // FIXME: you can use the bin width
      Ex.push_back (ex);
    }
  }

  // Put the fit results in arrays for TGraphErrors
  // ----------------------------------------------
  const int nn = means.size();
  double *x    = new double[nn];
  double *ym   = new double[nn];
  double *e    = new double[nn];
  double *eym  = new double[nn];
  double *yw   = new double[nn];
  double *eyw  = new double[nn];
  double *yc   = new double[nn];

  for (int j=0; j<nn; j++) {
    x[j]   = Xcenter[j];
    ym[j]  = means[j];
    eym[j] = Emeans[j];
    // yw[j]  = maxs[j];
    // eyw[j] = Emaxs[j];
    yw[j]  = sigmas[j];
    eyw[j] = Esigmas[j];
    yc[j]  = chi2s[j];
    e[j]   = Ex[j];
  }

  // Create TGraphErrors
  // -------------------
  TString name = histo->GetName();
  TGraphErrors *grM = new TGraphErrors (nn, x, ym, e, eym);
  grM->SetTitle (name+"_mean");
  grM->SetName (name+"_mean");
  TGraphErrors *grW = new TGraphErrors (nn, x, yw, e, eyw);
  grW->SetTitle (name+"_sigma");
  grW->SetName (name+"_sigma");
  TGraphErrors *grC = new TGraphErrors (nn, x, yc, e, e);
  grC->SetTitle (name+"_chi2");
  grC->SetName (name+"_chi2");

  // Cleanup
  // -------
  delete x;
  delete ym;
  delete eym;
  delete yw;
  delete eyw;
  delete yc;
  delete e;
  delete fitFcn;

  results.push_back (grM);
  results.push_back (grW);
  results.push_back (grC);
  return results;
}

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

Definition at line 476 of file MuScleFitUtils.cc.

References funct::cos(), funct::exp(), muMass, funct::sin(), mathSSE::sqrt(), and tmp.

Referenced by applyBias(), applyScale(), and applySmearing().

{
  double px = ptEtaPhiE[0]*cos(ptEtaPhiE[2]);
  double py = ptEtaPhiE[0]*sin(ptEtaPhiE[2]);
  double tmp = 2*atan(exp(-ptEtaPhiE[1]));
  double pz = ptEtaPhiE[0]*cos(tmp)/sin(tmp);
  double E  = sqrt(px*px+py*py+pz*pz+muMass*muMass);

  // lorentzVector corrMu(px,py,pz,E);
  // To fix memory leaks, this is to be substituted with
  // std::auto_ptr<lorentzVector> corrMu(new lorentzVector(px, py, pz, E));

  return lorentzVector(px,py,pz,E);
}

double MuScleFitUtils::invDimuonMass	(	const lorentzVector &	mu1,
		const lorentzVector &	mu2
	)		`[inline, static]`

Definition at line 493 of file MuScleFitUtils.cc.

Referenced by likelihood(), and minimizeLikelihood().

{
  return (mu1+mu2).mass();
}

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

Definition at line 2068 of file MuScleFitUtils.cc.

References gather_cfg::cout, debug, ExpressReco_HICollisions_FallBack::e, GL, runTheMatrix::np, P, Pi, ResGamma, ResMass, and x.

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

{
  // This routine computes the likelihood that a given measured mass "measMass" is
  // the result of resonance #ires if the resolution expected for the two muons is massResol
  // ---------------------------------------------------------------------------------------

  double P = 0.;

  // Return Lorentz value for zero resolution cases (like MC)
  // --------------------------------------------------------
  if (massResol==0.) {
    if (debug>9) std::cout <<  "Mass, gamma , mref, width, P: " << mass
                      << " " << ResGamma[ires] << " " << ResMass[ires]<< " " << massResol
                      << " : used Lorentz P-value" << std::endl;
    return (0.5*ResGamma[ires]/TMath::Pi())/((mass-ResMass[ires])*(mass-ResMass[ires])+
                                             .25*ResGamma[ires]*ResGamma[ires]);
  }

  // NB defined as below, P is not a "probability" but a likelihood that we observe
  // a dimuon mass "mass", given mRef, gamma, and massResol. It is what we need for the
  // fit which finds the best resolution parameters, though. A definition which is
  // more properly a probability is given below (in massProb2()), where the result
  // cannot be used to fit resolution parameters because the fit would always prefer
  // to set the res parameters to the minimum possible value (best resolution),
  // to have a probability close to one of observing any mass.
  // -------------------------------------------------------------------------------
  // NNBB the following two lines have been replaced with the six following them,
  // which provide an improvement of a factor 9 in speed of execution at a
  // negligible price in precision.
  // ----------------------------------------------------------------------------
  // GL->SetParameters(gamma,mRef,mass,massResol);
  // P = GL->Integral(mass-5*massResol, mass+5*massResol);

  Int_t np = 100;
  double * x = new double[np];
  double * w = new double[np];
  GL->SetParameters (ResGamma[ires], ResMass[ires], mass, massResol);
  GL->CalcGaussLegendreSamplingPoints (np, x, w, 0.1e-15);
  P = GL->IntegralFast (np, x, w, ResMass[ires]-10*ResGamma[ires], ResMass[ires]+10*ResGamma[ires]);
  delete[] x;
  delete[] w;

  // If we are too far away we set P to epsilon and forget about this event
  // ----------------------------------------------------------------------
  if (P<1.0e-12) {
    P = 1.0e-12;
    if (debug>9) std::cout << "Mass, gamma , mref, width, P: " << mass
                      << " " << ResGamma[ires] << " " << ResMass[ires] << " " << massResol
                      << ": used epsilon" << std::endl;
    return P;
  }

  if (debug>9) std::cout << "Mass, gamma , mref, width, P: " << mass
                    << " " << ResGamma[ires] << " " << ResMass[ires] << " " << massResol
                    << " " << P << std::endl;
  return P;
}

double MuScleFitUtils::massProb	(	const double &	mass,
		const double &	resEta,
		const double &	rapidity,
		const double &	massResol,
		const std::vector< double > &	parval,
		const bool	doUseBkgrWindow = `false`
	)		`[static]`

Definition at line 686 of file MuScleFitUtils.cc.

References ExpressReco_HICollisions_FallBack::id, massProb(), and L1TEmulatorMonitor_cff::p.

{
#ifdef USE_CALLGRIND
  CALLGRIND_START_INSTRUMENTATION;
#endif

  double * p = new double[(int)(parval.size())];
  // Replaced by auto_ptr, which handles delete at the end
  // Removed auto_ptr, check massResolution for an explanation.
  // std::auto_ptr<double> p(new double[(int)(parval.size())]);

  std::vector<double>::const_iterator it = parval.begin();
  int id = 0;
  for ( ; it!=parval.end(); ++it, ++id) {
    // (&*p)[id] = *it;
    p[id] = *it;
  }
  // p must be passed by value as below:
  double massProbability = massProb( mass, resEta, rapidity, massResol, p, doUseBkgrWindow );
  delete[] p;

#ifdef USE_CALLGRIND
  CALLGRIND_STOP_INSTRUMENTATION;
  CALLGRIND_DUMP_STATS;
#endif

  return massProbability;
}

double MuScleFitUtils::massProb	(	const double &	mass,
		const double &	resEta,
		const double &	rapidity,
		const double &	massResol,
		double *	parval,
		const bool	doUseBkgrWindow = `false`
	)		`[static]`

Definition at line 851 of file MuScleFitUtils.cc.

References BackgroundHandler::backgroundFunction(), backgroundHandler, backgroundProb_, checkMassWindow(), gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, GLNorm, GLValue, GLZNorm, GLZValue, i, j, loopCounter, minuitLoop_, MuonType, P, parBgr, CrossSectionHandler::parNum(), parResol, parScale, probability(), rapidityBinsForZ_, CrossSectionHandler::relativeCrossSections(), resfind, ResHalfWidth, BackgroundHandler::resMass(), ResMass, signalProb_, totalResNum, and BackgroundHandler::windowBorders().

                                                                                                                                                                            {

  // This routine computes the likelihood that a given measured mass "measMass" is
  // the result of a reference mass ResMass[] if the resolution
  // expected for the two muons is massResol.
  // This version includes two parameters (the last two in parval, by default)
  // to size up the background fraction and its relative normalization with respect
  // to the signal shape.
  //
  // We model the signal probability with a Lorentz L(M,H) of resonance mass M and natural width H
  // convoluted with a gaussian G(m,s) of measured mass m and expected mass resolution s,
  // by integrating over the intersection of the supports of L and G (which can be made to coincide with
  // 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:
  //
  //   GL(m,s) = Int(M-10H,M+10H) [ L(x-M,H) * G(x-m,s) ] dx
  //
  // The above convolution is computed numerically by an independent root macro, Probs.C, which outputs
  // the values in six 1001x1001 grids, one per resonance.
  //
  // NB THe following block of explanations for background models is outdated, see detailed
  // explanations where the code computes PB.
  // +++++++++++++++++++++++
  //   For the background, instead, we have two choices: a linear and an exponential model.
  //     * For the linear model, we choose a one-parameter form whereby the line is automatically normalized
  //       in the support [x1,x2] where we defined our "signal region", as follows:
  //
  //         B(x;b) = 1/(x2-x1) + {x - (x2+x1)/2} * b
  //
  //       Defined as above, B(x) is a line passing for the point of coordinates (x1+x2)/2, 1/(x2-x1),
  //       whose slope b has as support the interval ( -2/[(x1-x2)*(x1+x2)], 2/[(x1-x2)*(x1+x2)] )
  //       so that B(x) is always positive-definite in [x1,x2].
  //
  //     * For the exponential model, we define B(x;b) as
  //
  //         B(x;b) = b * { exp(-b*x) / [exp(-b*x1)-exp(-b*x2)] }
  //
  //       This one-parameter definition is automatically normalized to unity in [x1,x2], with a parameter
  //       b which has to be positive in order for the slope to be negative.
  //       Please note that this model is not useful in most circumstances; a more useful form would be one
  //       which included a linear component.
  // ++++++++++++++++++++++
  //
  // Once GL(m,s) and B(x;b) are defined, we introduce a further parameter a, such that we can have the
  // likelihood control the relative fraction of signal and background. We first normalize GL(m,s) for
  // any given s by taking the integral
  //
  //   Int(x1,x2) GL(m,s) dm = K_s
  //
  // We then define the probability as
  //
  //   P(m,s,a,b) = GL(m,s)/K_s * a  +  B(x,b) * (1-a)
  //
  // with a taking on values in the interval [0,1].
  // Defined as above, the probability is well-behaved, in the sense that it has a value between 0 and 1,
  // and the four parameters m,s,a,b fully control its shape.
  //
  // It is to be noted that the formulation above requires the computation of two rather time-consuming
  // integrals. The one defining GL(m,s) can be stored in a TH2D and loaded by the constructor from a
  // file suitably prepared, and this will save loads of computing time.
  // ----------------------------------------------------------------------------------------------------

  double P = 0.;
  int crossSectionParShift = parResol.size() + parScale.size();
  // Take the relative cross sections
  std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parval[crossSectionParShift]), resfind);
  //   for( unsigned int i=0; i<relativeCrossSections.size(); ++i ) {
  //     std::cout << "relativeCrossSections["<<i<<"] = " << relativeCrossSections[i] << std::endl;
  //     std::cout << "parval["<<crossSectionParShift+i<<"] = " << parval[crossSectionParShift+i] << std::endl;
  //   }

  // int bgrParShift = crossSectionParShift + parCrossSection.size();
  int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
  double Bgrp1 = 0.;
  //   double Bgrp2 = 0.;
  //   double Bgrp3 = 0.;

  // NB defined as below, P is a non-rigorous "probability" that we observe
  // a dimuon mass "mass", given ResMass[], gamma, and massResol. It is what we need for the
  // fit which finds the best resolution parameters, though. A definition which is
  // more properly a probability is given below (in massProb2()), where the result
  // cannot be used to fit resolution parameters because the fit would always prefer
  // to set the res parameters to the minimum possible value (best resolution),
  // to have a probability close to one of observing any mass.
  // -------------------------------------------------------------------------------

  // Determine what resonance(s) we have to deal with
  // NB for now we assume equal xs for each resonance
  // so we do not assign them different weights
  // ------------------------------------------------
  double PS[6] = {0.};
  double PB = 0.;
  double PStot[6] = {0.};

  // Should be removed because it is not used
  bool resConsidered[6] = {false};

  bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
  // bool useBackgroundWindow = (doBackgroundFit[loopCounter]);

  // First check the Z, which is divided in 24 rapidity bins
  // NB max value of Z rapidity to be considered is 2.4 here
  // -------------------------------------------------------

  // Do this only if we want to use the rapidity bins for the Z
  if( MuScleFitUtils::rapidityBinsForZ_ ) {
    // ATTENTION: cut on Z rapidity at 2.4 since we only have histograms up to that value
    // std::pair<double, double> windowFactors = backgroundHandler->windowFactors( useBackgroundWindow, 0 );
    std::pair<double, double> windowBorders = backgroundHandler->windowBorders( useBackgroundWindow, 0 );
    if( resfind[0]>0
        // && checkMassWindow( mass, 0,
        //                     backgroundHandler->resMass( useBackgroundWindow, 0 ),
        //                     windowFactors.first, windowFactors.second )
        && checkMassWindow( mass, windowBorders.first, windowBorders.second )
        // && fabs(rapidity)<2.4
        ) {

      int iY = (int)(fabs(rapidity)*10.);
      if( iY > 23 ) iY = 23;

      if (MuScleFitUtils::debug>1) std::cout << "massProb:resFound = 0, rapidity bin =" << iY << std::endl;

      // In this case the last value is the rapidity bin
      PS[0] = probability(mass, massResol, GLZValue, GLZNorm, 0, iY);

      if( PS[0] != PS[0] ) {
        std::cout << "ERROR: PS[0] = nan, setting it to 0" << std::endl;
        PS[0] = 0;
      }

      std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
                                                                                   &(parval[bgrParShift]), MuScleFitUtils::totalResNum, 0,
                                                                                   resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
      Bgrp1 = bgrResult.first;
      // When fitting the background we have only one Bgrp1
      // When not fitting the background we have many only in a superposition region and this case is treated
      // separately after this loop
      PB = bgrResult.second;
      if( PB != PB ) PB = 0;
      PStot[0] = (1-Bgrp1)*PS[0] + Bgrp1*PB;

      // PStot[0] *= crossSectionHandler->crossSection(0);
      // PStot[0] *= parval[crossSectionParShift];
      PStot[0] *= relativeCrossSections[0];
      // std::cout << "PStot["<<0<<"] = " << "(1-"<<Bgrp1<<")*"<<PS[0]<<" + "<<Bgrp1<<"*"<<PB<<" = " << PStot[0] << std::endl;
    }
    else {
      if( debug > 0 ) {
        std::cout << "Mass = " << mass << " outside range with rapidity = " << rapidity << std::endl;
        std::cout << "with resMass = " << backgroundHandler->resMass( useBackgroundWindow, 0 ) << " and left border = " << windowBorders.first << " right border = " << windowBorders.second << std::endl;
      }
    }
  }
  // Next check the other resonances
  // -------------------------------
  int firstRes = 1;
  if( !MuScleFitUtils::rapidityBinsForZ_ ) firstRes = 0;
  for( int ires=firstRes; ires<6; ++ires ) {
    if( resfind[ires] > 0 ) {
      // First is left, second is right (returns (1,1) in the case of resonances, it could be improved avoiding the call in this case)
      // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
      std::pair<double, double> windowBorder = backgroundHandler->windowBorders( useBackgroundWindow, ires );
      if( checkMassWindow(mass, windowBorder.first, windowBorder.second) ) {
        if (MuScleFitUtils::debug>1) std::cout << "massProb:resFound = " << ires << std::endl;

        // In this case the rapidity value is instead the resonance index again.
        PS[ires] = probability(mass, massResol, GLValue, GLNorm, ires, ires);

        std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
                                                                                     &(parval[bgrParShift]), MuScleFitUtils::totalResNum, ires,
                                                                                     resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
        Bgrp1 = bgrResult.first;
        PB = bgrResult.second;

        if( PB != PB ) PB = 0;
        PStot[ires] = (1-Bgrp1)*PS[ires] + Bgrp1*PB;
        if( MuScleFitUtils::debug>0 ) std::cout << "PStot["<<ires<<"] = " << "(1-"<<Bgrp1<<")*"<<PS[ires]<<" + "<<Bgrp1<<"*"<<PB<<" = " << PStot[ires] << std::endl;

        PStot[ires] *= relativeCrossSections[ires];
      }
    }
  }

  for( int i=0; i<6; ++i ) {
    P += PStot[i];
  }

  if( MuScleFitUtils::signalProb_ != 0 && MuScleFitUtils::backgroundProb_ != 0 ) {
    double PStotTemp = 0.;
    for( int i=0; i<6; ++i ) {
      PStotTemp += PS[i]*relativeCrossSections[i];
    }
    if( PStotTemp != PStotTemp ) {
      std::cout << "ERROR: PStotTemp = nan!!!!!!!!!" << std::endl;
      int parnumber = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size());
      for( int i=0; i<6; ++i ) {
        std::cout << "PS[i] = " << PS[i] << std::endl;
        if( PS[i] != PS[i] ) {
          std::cout << "mass = " << mass << std::endl;
          std::cout << "massResol = " << massResol << std::endl;
          for( int j=0; j<parnumber; ++j ) {
            std::cout << "parval["<<j<<"] = " << parval[j] << std::endl;
          }
        }
      }
    }
    if( PStotTemp == PStotTemp ) {
      MuScleFitUtils::signalProb_->SetBinContent(MuScleFitUtils::minuitLoop_, MuScleFitUtils::signalProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PStotTemp);
    }
    if (debug>0) std::cout << "mass = " << mass << ", P = " << P << ", PStot = " << PStotTemp << ", PB = " << PB << ", bgrp1 = " << Bgrp1 << std::endl;

    MuScleFitUtils::backgroundProb_->SetBinContent(MuScleFitUtils::minuitLoop_, MuScleFitUtils::backgroundProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PB);
  }
  return P;
}

double MuScleFitUtils::massResolution	(	const lorentzVector &	mu1,
		const lorentzVector &	mu2,
		double *	parval
	)		`[static]`

Definition at line 528 of file MuScleFitUtils.cc.

References 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, funct::exp(), LogDebug, massResolComponents, mMu2, funct::pow(), resfind, ResHalfWidth, ResMass, ResMaxSigma, resolutionFunction, resolutionFunctionBase< T >::sigmaCotgTh(), resolutionFunctionBase< T >::sigmaPhi(), resolutionFunctionBase< T >::sigmaPt(), funct::sin(), and mathSSE::sqrt().

{
  // 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.
  // --------------------------------------------------------------------------
  double mass   = (mu1+mu2).mass();
  double pt1    = mu1.Pt();
  double phi1   = mu1.Phi();
  double eta1   = mu1.Eta();
  double theta1 = 2*atan(exp(-eta1));
  // double cotgTheta1 = cos(theta1)/sin(theta1);
  double pt2    = mu2.Pt();
  double phi2   = mu2.Phi();
  double eta2   = mu2.Eta();
  double theta2 = 2*atan(exp(-eta2));
  // double cotgTheta2 = cos(theta2)/sin(theta2);

  // double mass_check = sqrt(2*mMu2+2*sqrt(std::pow(pt1/sin(theta1),2)+mMu2)*sqrt(std::pow(pt2/sin(theta2),2)+mMu2)-
  //                       2*pt1*pt2*(cos(phi1-phi2)+1/(tan(theta1)*tan(theta2))));

  // ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
  // -----------------------------------------------------------------------------------------------------------
  double dmdpt1  = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))-
                    pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
  double dmdpt2  = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))-
                    pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
  double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
  double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
  // double dmdtheta1 = (-std::pow(pt1/sin(theta1),2)/tan(theta1)*
  //                    sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))+
  //                    2*pt1*pt2/(tan(theta2)*std::pow(sin(theta1),2)))/mass;
  // double dmdtheta2 = (-std::pow(pt2/sin(theta2),2)/tan(theta2)*
  //                    sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))+
  //                    2*pt2*pt1/(tan(theta1)*std::pow(sin(theta2),2)))/mass;
  double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
                       sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) -
                       pt1*pt2*cos(theta2)/sin(theta2))/mass;
  double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
                       sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) -
                       pt2*pt1*cos(theta1)/sin(theta1))/mass;

  if( debugMassResol_ ) {
    massResolComponents.dmdpt1 = dmdpt1;
    massResolComponents.dmdpt2 = dmdpt2;
    massResolComponents.dmdphi1 = dmdphi1;
    massResolComponents.dmdphi2 = dmdphi2;
    massResolComponents.dmdcotgth1 = dmdcotgth1;
    massResolComponents.dmdcotgth2 = dmdcotgth2;
  }

  // Resolution parameters:
  // ----------------------
  double sigma_pt1 = resolutionFunction->sigmaPt( pt1,eta1,parval );
  double sigma_pt2 = resolutionFunction->sigmaPt( pt2,eta2,parval );
  double sigma_phi1 = resolutionFunction->sigmaPhi( pt1,eta1,parval );
  double sigma_phi2 = resolutionFunction->sigmaPhi( pt2,eta2,parval );
  double sigma_cotgth1 = resolutionFunction->sigmaCotgTh( pt1,eta1,parval );
  double sigma_cotgth2 = resolutionFunction->sigmaCotgTh( pt2,eta2,parval );
  double cov_pt1pt2 = resolutionFunction->covPt1Pt2( pt1, eta1, pt2, eta2, parval );

  // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
  // multiply it by pt.
  double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2)+
                         std::pow(dmdphi1*sigma_phi1,2)+std::pow(dmdphi2*sigma_phi2,2)+
                         std::pow(dmdcotgth1*sigma_cotgth1,2)+std::pow(dmdcotgth2*sigma_cotgth2,2)+
                         2*dmdpt1*dmdpt2*cov_pt1pt2);

  // if( sigma_cotgth1 < 0 || sigma_cotgth2 < 0 ) {
  //   std::cout << "WARNING: sigma_cotgth1 = " << sigma_cotgth1 << std::endl;
  //   std::cout << "WARNING: sigma_cotgth2 = " << sigma_cotgth2 << std::endl;
  //   std::cout << "mass_res = " << mass_res << std::endl;
  // }

  // double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2));

  if (debug>19) {
    std::cout << "  Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
         << " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
    std::cout << " P[0]="
         << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << std::endl;
    std::cout << "  Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1="
         << sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
    std::cout << "  Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1="
         << sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
    std::cout << "  Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2
         << " sigma_cotgth1="
         << sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
    std::cout << "  Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2
         << " : " << mass << " +- " << mass_res << std::endl;
  }

  // Debug std::cout
  // ----------
  bool didit = false;
  for (int ires=0; ires<6; ires++) {
    if (!didit && resfind[ires]>0 && fabs(mass-ResMass[ires])<ResHalfWidth[ires]) {
      if (mass_res>ResMaxSigma[ires] && counter_resprob<100) {
        counter_resprob++;
        LogDebug("MuScleFitUtils") << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
        // std::cout << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
//      std::cout << "---------------------------" << std::endl;
//      std::cout << "  Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
//           << " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
//      if (ResolFitType==1)
//        std::cout << " P[0]="
//             << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3]
//             << " P[4]=" << parval[4] << " P[5]=" << parval[5] << std::endl;
//      if (ResolFitType==2)
//        std::cout << " P[0]="
//             << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << std::endl;
//      std::cout << "  Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1="
//           << sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
//      std::cout << "  Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1="
//           << sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
//      std::cout << "  Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2
//           << " sigma_cotgth1="
//           << sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
//      std::cout << "  Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2
//           << " : " << mass << " +- " << mass_res << std::endl;
//      std::cout << "---------------------------" << std::endl;
        didit = true;
      }
    }
  }

//   if( mass_res != mass_res ) {
//     std::cout << "MASS_RESOL_NAN PROBLEM:" << std::endl;
//     std::cout << "---------------------------" << std::endl;
//     std::cout << "  Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
//          << " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
//     std::cout << " P[0]=" << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << " P[4]=" << parval[4] << std::endl;
//     std::cout << "  Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1=" << sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
//     std::cout << "  Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1=" << sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
//     std::cout << "  Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2 << " sigma_cotgth1=" << sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
//     std::cout << "  Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2 << " : " << mass << " +- " << mass_res << std::endl;
//     std::cout << "---------------------------" << std::endl;
//   }

  return mass_res;
}

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

Definition at line 501 of file MuScleFitUtils.cc.

References ExpressReco_HICollisions_FallBack::id, massResolution(), and L1TEmulatorMonitor_cff::p.

{
  // double * p = new double[(int)(parval.size())];
  // Replaced by auto_ptr, which handles delete at the end
  // --------- //
  // ATTENTION //
  // --------- //
  // auto_ptr calls delete, not delete[] and thus it must
  // not be used with arrays. There are alternatives see
  // e.g.: http://www.gotw.ca/gotw/042.htm. The best
  // alternative seems to be to switch to vector though.
  // std::auto_ptr<double> p(new double[(int)(parval.size())]);

  double * p = new double[(int)(parval.size())];
  std::vector<double>::const_iterator it = parval.begin();
  int id = 0;
  for ( ; it!=parval.end(); ++it, ++id) {
    // (&*p)[id] = *it;
    p[id] = *it;
  }
  double massRes = massResolution (mu1, mu2, p);
  delete[] p;
  return massRes;
}

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

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

static double MuScleFitUtils::massResolution	(	const lorentzVector &	mu1,
		const lorentzVector &	mu2,
		std::auto_ptr< double >	parval
	)		`[static]`

void MuScleFitUtils::minimizeLikelihood ( ) [static]

Definition at line 1112 of file MuScleFitUtils.cc.

References backgroundHandler, backgroundProb_, BgrFitType, svgfig::canvas(), CastorDataFrameFilter_impl::check(), checkMassWindow(), computeMinosErrors_, gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, doCrossSectionFit, doResolFit, doScaleFit, duringMinos_, FitStrategy, i, invDimuonMass(), likelihood(), likelihoodInLoop_, loopCounter, massWindowHalfWidth, minimumShapePlots_, minuitLoop_, MuonType, AlCaRecoCosmics_cfg::name, normalizationChanged_, parBgr, parBgrFix, parBgrOrder, parCrossSection, parCrossSectionOrder, parfix, scaleFunctionBase< T >::parNum(), CrossSectionHandler::parNum(), resolutionFunctionBase< T >::parNum(), parorder, parResol, parResolFix, parResolOrder, parScale, parScaleFix, parScaleOrder, parvalue, ReducedSavedPair, BackgroundHandler::regionsParNum(), CrossSectionHandler::relativeCrossSections(), CrossSectionHandler::releaseParameters(), BackgroundHandler::rescale(), scaleFunctionBase< T >::resetParameters(), resfind, BackgroundHandler::resMass(), ResMass, ResolFitType, rminPtr_, SavedPair, scaleFitNotDone_, ScaleFitType, scaleFunction, CrossSectionHandler::setParameters(), BackgroundHandler::setParameters(), resolutionFunctionBase< T >::setParameters(), scaleFunctionBase< T >::setParameters(), signalProb_, startWithSimplex_, cmsPerfCommons::Step, BackgroundHandler::unlockParameter(), and BackgroundHandler::windowBorders().

Referenced by MuScleFit::endOfFastLoop().

{
  // Output file with fit parameters resulting from minimization
  // -----------------------------------------------------------
  ofstream FitParametersFile;
  FitParametersFile.open ("FitParameters.txt", std::ios::app);
  FitParametersFile << "Fitting with resolution, scale, bgr function # "
                    << ResolFitType << " " << ScaleFitType << " " << BgrFitType
                    << " - Iteration " << loopCounter << std::endl;

  // Fill parvalue and other vectors needed for the fitting
  // ------------------------------------------------------





  // ----- //
  // FIXME //
  // ----- //
  // this was changed to verify the possibility that fixed parameters influence the errors.
  // It must be 0 otherwise the parameters for resonances will not be passed by minuit (will be always 0).
  // Should be removed.
  int parForResonanceWindows = 0;
  // int parnumber = (int)(parResol.size()+parScale.size()+parCrossSection.size()+parBgr.size() - parForResonanceWindows);
  int parnumber = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size() - parForResonanceWindows);






  int parnumberAll = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size());

  // parvalue is a std::vector<std::vector<double> > storing all the parameters from all the loops
  parvalue.push_back(parResol);
  std::vector<double> *tmpVec = &(parvalue.back());

  // If this is not the first loop we want to start from neutral values
  // Otherwise the scale will start with values correcting again a bias
  // that is already corrected.
  if( scaleFitNotDone_ ) {
    tmpVec->insert( tmpVec->end(), parScale.begin(), parScale.end() );
    std::cout << "scaleFitNotDone: tmpVec->size() = " << tmpVec->size() << std::endl;
  }
  else {
    scaleFunction->resetParameters(tmpVec);
    std::cout << "scaleFitDone: tmpVec->size() = " << tmpVec->size() << std::endl;
  }
  tmpVec->insert( tmpVec->end(), parCrossSection.begin(), parCrossSection.end() );
  tmpVec->insert( tmpVec->end(), parBgr.begin(), parBgr.end() );
  int i = 0;
  std::vector<double>::const_iterator it = tmpVec->begin();
  for( ; it != tmpVec->end(); ++it, ++i ) {
    std::cout << "tmpVec["<<i<<"] = " << *it << std::endl;
  }

  // Empty vector of size = number of cross section fitted parameters. Note that the cross section
  // fit works in a different way than the others and it uses ratios of the paramters passed via cfg.
  // We use this empty vector for compatibility with the rest of the structure.
  std::vector<int> crossSectionParNumSizeVec( MuScleFitUtils::crossSectionHandler->parNum(), 0 );

  std::vector<int> parfix(parResolFix);
  parfix.insert( parfix.end(), parScaleFix.begin(), parScaleFix.end() );
  parfix.insert( parfix.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end() );
  parfix.insert( parfix.end(), parBgrFix.begin(), parBgrFix.end() );

  std::vector<int> parorder(parResolOrder);
  parorder.insert( parorder.end(), parScaleOrder.begin(), parScaleOrder.end() );
  parorder.insert( parorder.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end() );
  parorder.insert( parorder.end(), parBgrOrder.begin(), parBgrOrder.end() );

  // This is filled later
  std::vector<double> parerr(3*parnumberAll,0.);

  if (debug>19) {
    std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters before likelihood " << std::endl;
    for (unsigned int i=0; i<(unsigned int)parnumberAll; i++) {
      std::cout << "  Par # " << i << " = " << parvalue[loopCounter][i] << " : free = " << parfix[i] << "; order = "
           << parorder[i] << std::endl;
    }
  }

//   // Background rescaling from regions to resonances
//   // -----------------------------------------------
//   // If we are in a loop > 0 and we are not fitting the background, but we have fitted it in the previous iteration
//   if( loopCounter > 0 && !(doBackgroundFit[loopCounter]) && doBackgroundFit[loopCounter-1] ) {
//     // This rescales from regions to resonances
//     int localMuonType = MuonType;
//     if( MuonType > 2 ) localMuonType = 2;
//     backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
//                                 MuScleFitUtils::SavedPair);
//   }

  // Init Minuit
  // -----------
  TMinuit rmin (parnumber);
  rminPtr_ = &rmin;
  rmin.SetFCN (likelihood);     // Unbinned likelihood
  // Standard initialization of minuit parameters:
  // sets input to be $stdin, output to be $stdout
  // and saving to a file.
  rmin.mninit (5, 6, 7);
  int ierror = 0;
  int istat;
  double arglis[4];
  arglis[0] = FitStrategy;      // Strategy 1 or 2
  // 1 standard
  // 2 try to improve minimum (slower)
  rmin.mnexcm ("SET STR", arglis, 1, ierror);

  arglis[0] = 10001;
  // Set the random seed for the generator used in SEEk to a fixed value for reproducibility
  rmin.mnexcm("SET RAN", arglis, 1, ierror);

  // Set fit parameters
  // ------------------
  double * Start = new double[parnumberAll];
  double * Step  = new double[parnumberAll];
  double * Mini  = new double[parnumberAll];
  double * Maxi  = new double[parnumberAll];
  int * ind = new int[parnumberAll]; // Order of release of parameters
  TString * parname = new TString[parnumberAll];

  MuScleFitUtils::resolutionFunctionForVec->setParameters( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, MuonType );

  // Take the number of parameters in the resolutionFunction and displace the arrays passed to the scaleFunction
  int resParNum = MuScleFitUtils::resolutionFunctionForVec->parNum();

  MuScleFitUtils::scaleFunctionForVec->setParameters( &(Start[resParNum]), &(Step[resParNum]), &(Mini[resParNum]), &(Maxi[resParNum]),
                                                      &(ind[resParNum]), &(parname[resParNum]), parScale, parScaleOrder, MuonType );

  // Initialize cross section parameters
  int crossSectionParShift = resParNum + MuScleFitUtils::scaleFunctionForVec->parNum();
  MuScleFitUtils::crossSectionHandler->setParameters( &(Start[crossSectionParShift]), &(Step[crossSectionParShift]), &(Mini[crossSectionParShift]),
                                                      &(Maxi[crossSectionParShift]), &(ind[crossSectionParShift]), &(parname[crossSectionParShift]),
                                                      parCrossSection, parCrossSectionOrder, resfind );

  // Initialize background parameters
  int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
  MuScleFitUtils::backgroundHandler->setParameters( &(Start[bgrParShift]), &(Step[bgrParShift]), &(Mini[bgrParShift]), &(Maxi[bgrParShift]),
                                                    &(ind[bgrParShift]), &(parname[bgrParShift]), parBgr, parBgrOrder, MuonType );

  for( int ipar=0; ipar<parnumber; ++ipar ) {
    std::cout << "parname["<<ipar<<"] = " << parname[ipar] << std::endl;
    std::cout << "Start["<<ipar<<"] = " << Start[ipar] << std::endl;
    std::cout << "Step["<<ipar<<"] = " << Step[ipar] << std::endl;
    std::cout << "Mini["<<ipar<<"] = " << Mini[ipar] << std::endl;
    std::cout << "Maxi["<<ipar<<"] = " << Maxi[ipar] << std::endl;


    rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], Mini[ipar], Maxi[ipar], ierror );


    // Testing without limits
    // rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], 0, 0, ierror );


  }

  // Do minimization
  // ---------------
  if (debug>19)
    std::cout << "[MuScleFitUtils-minimizeLikelihood]: Starting minimization" << std::endl;
  double fmin;
  double fdem;
  double errdef;
  int npari;
  int nparx;
  rmin.mnexcm ("call fcn", arglis, 1, ierror);

  // First, fix all parameters
  // -------------------------
  if (debug>19)
    std::cout << "[MuScleFitUtils-minimizeLikelihood]: First fix all parameters ...";
  for (int ipar=0; ipar<parnumber; ipar++) {
    rmin.FixParameter (ipar);
  }

  // Then release them in the specified order and refit
  // --------------------------------------------------
  if (debug>19) std::cout << " Then release them in order..." << std::endl;

  TString name;
  double pval;
  double pmin;
  double pmax;
  double errp;
  double errl;
  double errh;
  int ivar;
  double erro;
  double cglo;
  int n_times = 0;
  // n_times = number of loops required to unlock all parameters.

  if (debug>19) std::cout << "Before scale parNum" << std::endl;
  int scaleParNum = scaleFunction->parNum();
  if (debug>19) std::cout << "After scale parNum" << std::endl;
  // edm::LogInfo("minimizeLikelihood") << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
  // edm::LogInfo("minimizeLikelihood") << "number of parameters for resolutionFunction = " << resParNum << std::endl;
  // edm::LogInfo("minimizeLikelihood") << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
  // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
  std::cout << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
  std::cout << "number of parameters for resolutionFunction = " << resParNum << std::endl;
  std::cout << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
  std::cout << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
  // edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << backgroundFunction->parNum() << std::endl;

  for (int i=0; i<parnumber; i++) {
    // NB ind[] has been set as parorder[] previously
    if (n_times<ind[i]) {
      edm::LogInfo("minimizeLikelihood") << "n_times = " << n_times << ", ind["<<i<<"] = " << ind[i] << ", scaleParNum = " << scaleParNum << ", doScaleFit["<<loopCounter<<"] = " << doScaleFit[loopCounter] << std::endl;
      // Set the n_times only if we will do the fit
      if ( i<resParNum ) {
        if( doResolFit[loopCounter] ) n_times = ind[i];
      }
      else if( i<resParNum+scaleParNum ) {
        if( doScaleFit[loopCounter] ) n_times = ind[i];
      }
      else if( doBackgroundFit[loopCounter] ) n_times = ind[i];
    }
  }
  for (int iorder=0; iorder<n_times+1; iorder++) { // Repeat fit n_times times
    std::cout << "Starting minimization " << iorder << " of " << n_times << std::endl;

    bool somethingtodo = false;

    // Use parameters from cfg to select which fit to do
    // -------------------------------------------------
    if( doResolFit[loopCounter] ) {
      // Release resolution parameters and fit them
      // ------------------------------------------
      for( unsigned int ipar=0; ipar<parResol.size(); ++ipar ) {
        if( parfix[ipar]==0 && ind[ipar]==iorder ) {
          rmin.Release( ipar );
          somethingtodo = true;
        }
      }
    }
    if( doScaleFit[loopCounter] ) {
      // Release scale parameters and fit them
      // -------------------------------------
      for( unsigned int ipar=parResol.size(); ipar<parResol.size()+parScale.size(); ++ipar ) {
        if( parfix[ipar]==0 && ind[ipar]==iorder ) { // parfix=0 means parameter is free
          rmin.Release( ipar );
          somethingtodo = true;
        }
      }
      scaleFitNotDone_ = false;
    }
    unsigned int crossSectionParShift = parResol.size()+parScale.size();
    if( doCrossSectionFit[loopCounter] ) {
      // Release cross section parameters and fit them
      // ---------------------------------------------
      // Note that only cross sections of resonances that are being fitted are released
      bool doCrossSection = crossSectionHandler->releaseParameters( rmin, resfind, parfix, ind, iorder, crossSectionParShift );
      if( doCrossSection ) somethingtodo = true;
    }
    if( doBackgroundFit[loopCounter] ) {
      // Release background parameters and fit them
      // ------------------------------------------
      // for( int ipar=parResol.size()+parScale.size(); ipar<parnumber; ++ipar ) {
      // Free only the parameters for the regions, as the resonances intervals are never used to fit the background
      unsigned int bgrParShift = crossSectionParShift+crossSectionHandler->parNum();
      for( unsigned int ipar = bgrParShift; ipar < bgrParShift+backgroundHandler->regionsParNum(); ++ipar ) {
        // Release only those parameters for the resonances we are fitting
        if( parfix[ipar]==0 && ind[ipar]==iorder && backgroundHandler->unlockParameter(resfind, ipar - bgrParShift) ) {
          rmin.Release( ipar );
          somethingtodo = true;
        }
      }
    }

    // OK, now do minimization if some parameter has been released
    // -----------------------------------------------------------
    if( somethingtodo ) {
// #ifdef DEBUG

      std::stringstream fileNum;
      fileNum << loopCounter;

      minuitLoop_ = 0;
      char name[50];
      sprintf(name, "likelihoodInLoop_%d_%d", loopCounter, iorder);
      TH1D * tempLikelihoodInLoop = new TH1D(name, "likelihood value in minuit loop", 10000, 0, 10000);
      likelihoodInLoop_ = tempLikelihoodInLoop;
      char signalProbName[50];
      sprintf(signalProbName, "signalProb_%d_%d", loopCounter, iorder);
      TH1D * tempSignalProb = new TH1D(signalProbName, "signal probability", 10000, 0, 10000);
      signalProb_ = tempSignalProb;
      char backgroundProbName[50];
      sprintf(backgroundProbName, "backgroundProb_%d_%d", loopCounter, iorder);
      TH1D * tempBackgroundProb = new TH1D(backgroundProbName, "background probability", 10000, 0, 10000);
      backgroundProb_ = tempBackgroundProb;
// #endif



      // Before we start the minimization we create a list of events with only the events inside a smaller
      // window than the one in which the probability is != 0. We will compute the probability for all those
      // events and hopefully the margin will avoid them to get a probability = 0 (which causes a discontinuity
      // in the likelihood function). The width of this smaller window must be optimized, but we can start using
      // an 90% of the normalization window.
      double protectionFactor = 0.9;

      MuScleFitUtils::ReducedSavedPair.clear();
      for( unsigned int nev=0; nev<MuScleFitUtils::SavedPair.size(); ++nev ) {
        const lorentzVector * recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
        const lorentzVector * recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
        double mass = MuScleFitUtils::invDimuonMass( *recMu1, *recMu2 );
        // Test all resonances to see if the mass is inside at least one of the windows
        bool check = false;
        for( int ires = 0; ires < 6; ++ires ) {
          // std::pair<double, double> windowFactor = backgroundHandler->windowFactors( doBackgroundFit[loopCounter], ires );
          std::pair<double, double> windowBorder = backgroundHandler->windowBorders( doBackgroundFit[loopCounter], ires );
          // if( resfind[ires] && checkMassWindow( mass, ires, backgroundHandler->resMass( doBackgroundFit[loopCounter], ires ),
          //                                       0.9*windowFactor.first, 0.9*windowFactor.second ) ) {
          double resMassValue = backgroundHandler->resMass( doBackgroundFit[loopCounter], ires );
          double windowBorderLeft = resMassValue - protectionFactor*(resMassValue - windowBorder.first);
          double windowBorderRight = resMassValue + protectionFactor*(windowBorder.second - resMassValue);
          if( resfind[ires] && checkMassWindow( mass, windowBorderLeft, windowBorderRight ) ) {
            check = true;
          }
        }
        if( check ) {
          MuScleFitUtils::ReducedSavedPair.push_back(std::make_pair(*recMu1, *recMu2));
        }
      }
      std::cout << "Fitting with " << MuScleFitUtils::ReducedSavedPair.size() << " events" << std::endl;


      // rmin.SetMaxIterations(500*parnumber);

      MuScleFitUtils::normalizationChanged_ = 0;

      // Maximum number of iterations
      arglis[0] = 5000;

      // Run simplex first to get an initial estimate of the minimum
      if( startWithSimplex_ ) {
        rmin.mnexcm( "simplex", arglis, 0, ierror );
      }

      rmin.mnexcm( "mini", arglis, 0, ierror );


// #ifdef DEBUG
      likelihoodInLoop_->Write();
      signalProb_->Write();
      backgroundProb_->Write();
      delete tempLikelihoodInLoop;
      delete tempSignalProb;
      delete tempBackgroundProb;
      likelihoodInLoop_ = 0;
      signalProb_ = 0;
      backgroundProb_ = 0;
// #endif


      // Compute again the error matrix
      rmin.mnexcm( "hesse", arglis, 0, ierror );

      // Peform minos error analysis.
      if( computeMinosErrors_ ) {
        duringMinos_ = true;
        rmin.mnexcm( "minos", arglis, 0, ierror );
        duringMinos_ = false;
      }

      if( normalizationChanged_ > 1 ) {
        std::cout << "WARNING: normalization changed during fit meaning that events exited from the mass window. This causes a discontinuity in the likelihood function. Please check the scan of the likelihood as a function of the parameters to see if there are discontinuities around the minimum." << std::endl;
      }
    }

    // bool notWritten = true;
    for (int ipar=0; ipar<parnumber; ipar++) {

      rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
      // Save parameters in parvalue[] vector
      // ------------------------------------
      if (ierror!=0 && debug>0) {
        std::cout << "[MuScleFitUtils-minimizeLikelihood]: ierror!=0, bogus pars" << std::endl;
      }
      //     for (int ipar=0; ipar<parnumber; ipar++) {
      //       rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
      parvalue[loopCounter][ipar] = pval;
      //     }


      // int ilax2 = 0;
      // Double_t val2pl, val2mi;
      // rmin.mnmnot (ipar+1, ilax2, val2pl, val2mi);
      rmin.mnerrs (ipar, errh, errl, errp, cglo);


      // Set error on params
      // -------------------
      if (errp!=0) {
        parerr[3*ipar] = errp;
      } else {
        parerr[3*ipar] = (((errh)>(fabs(errl)))?(errh):(fabs(errl)));
      }
      parerr[3*ipar+1] = errl;
      parerr[3*ipar+2] = errh;

      if( ipar == 0 ) {
        FitParametersFile << " Resolution fit parameters:" << std::endl;
      }
      if( ipar == int(parResol.size()) ) {
        FitParametersFile << " Scale fit parameters:" << std::endl;
      }
      if( ipar == int(parResol.size()+parScale.size()) ) {
        FitParametersFile << " Cross section fit parameters:" << std::endl;
      }
      if( ipar == int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) ) {
        FitParametersFile << " Background fit parameters:" << std::endl;
      }
//       if( ipar >= int(parResol.size()+parScale.size()) && ipar < int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) && notWritted ) {

//         std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]));
//         std::vector<double>::const_iterator it = relativeCrossSections.begin();
//         for( ; it != relativeCrossSections.end(); ++it ) {
//           FitParametersFile << "  Results of the fit: parameter " << ipar << " has value "
//                             << *it << "+-" << 0
//                             << " + " << 0 << " - " << 0
//                             << " /t/t (" << 0 << ")" << std::endl;
//         }

//         notWritten = false;
//       }
//       else {
        FitParametersFile << "  Results of the fit: parameter " << ipar << " has value "
                          << pval << "+-" << parerr[3*ipar]
                          << " + " << parerr[3*ipar+1] << " - " << parerr[3*ipar+2]
                          // << " \t\t (" << parname[ipar] << ")"
                          << std::endl;
//       }
    }
    rmin.mnstat (fmin, fdem, errdef, npari, nparx, istat); // NNBB Commented for a check!
    FitParametersFile << std::endl;

    if( minimumShapePlots_ ) {
      // Create plots of the minimum vs parameters
      // -----------------------------------------
      // Keep this after the parameters filling because it recomputes the values and it can compromise the fit results.
      if( somethingtodo ) {
        std::stringstream iorderString;
        iorderString << iorder;
        std::stringstream iLoopString;
        iLoopString << loopCounter;

        for (int ipar=0; ipar<parnumber; ipar++) {
          if( parfix[ipar] == 1 ) continue;
          std::cout << "plotting parameter = " << ipar+1 << std::endl;
          std::stringstream iparString;
          iparString << ipar+1;
          rmin.mncomd( ("scan "+iparString.str()).c_str(), ierror );
          if( ierror == 0 ) {
            TCanvas * canvas = new TCanvas(("likelihoodCanvas_loop_"+iLoopString.str()+"_oder_"+iorderString.str()+"_par_"+iparString.str()).c_str(), ("likelihood_"+iparString.str()).c_str(), 1000, 800);
            canvas->cd();
            // arglis[0] = ipar;
            // rmin.mnexcm( "sca", arglis, 0, ierror );
            TGraph * graph = (TGraph*)rmin.GetPlot();
            graph->Draw("AP");
            // graph->SetTitle(("parvalue["+iparString.str()+"]").c_str());
            graph->SetTitle(parname[ipar]);
            // graph->Write();

            canvas->Write();
          }
        }

        //       // Draw contours of the fit
        //       TCanvas * canvas = new TCanvas(("contourCanvas_oder_"+iorderString.str()).c_str(), "contour", 1000, 800);
        //       canvas->cd();
        //       TGraph * contourGraph = (TGraph*)rmin.Contour(4, 2, 4);
        //       if( (rmin.GetStatus() == 0) || (rmin.GetStatus() >= 3) ) {
        //         contourGraph->Draw("AP");
        //       }
        //       else {
        //         std::cout << "Contour graph error: status = " << rmin.GetStatus() << std::endl;
        //       }
        //       canvas->Write();
      }
    }

  } // end loop on iorder
  FitParametersFile.close();

  std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters after likelihood " << std::endl;
  for (unsigned int ipar=0; ipar<(unsigned int)parnumber; ipar++) {
    std::cout << ipar << " " << parvalue[loopCounter][ipar] << " : free = "
         << parfix[ipar] << "; order = " << parorder[ipar] << std::endl;
  }

  // Put back parvalue into parResol, parScale, parCrossSection, parBgr
  // ------------------------------------------------------------------
  for( int i=0; i<(int)(parResol.size()); ++i ) {
    parResol[i] = parvalue[loopCounter][i];
  }
  for( int i=0; i<(int)(parScale.size()); ++i ) {
    parScale[i] = parvalue[loopCounter][i+parResol.size()];
  }
  parCrossSection = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]), resfind);
  for( unsigned int i=0; i<parCrossSection.size(); ++i ) {
  //   parCrossSection[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()];
    std::cout << "relative cross section["<<i<<"] = " << parCrossSection[i] << std::endl;
  }
  // Save only the fitted background parameters
  for( unsigned int i = 0; i<(parBgr.size() - parForResonanceWindows); ++i ) {
    parBgr[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()+crossSectionHandler->parNum()];
  }

  // Background rescaling from regions to resonances
  // -----------------------------------------------
  // Only if we fitted the background
  if( doBackgroundFit[loopCounter] ) {
    // This rescales from regions to resonances
    int localMuonType = MuonType;
    if( MuonType > 2 ) localMuonType = 2;
    backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
                                MuScleFitUtils::ReducedSavedPair );
  }

  // Delete the arrays used to set some parameters
  delete[] Start;
  delete[] Step;
  delete[] Mini;
  delete[] Maxi;
  delete[] ind;
  delete[] parname;
}

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

References counter_resprob, gather_cfg::cout, debug, LogDebug, nbins, ResHalfWidth, ResMass, ResMaxSigma, and ResMinMass.

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

{
  if( iRes == 0 && iY > 23 ) {
    std::cout << "WARNING: rapidity bin selected = " << iY << " but there are only histograms for the first 24 bins" << std::endl;
  }

  double PS = 0.;
  bool insideProbMassWindow = true;
  // Interpolate the four values of GLZValue[] in the
  // grid square within which the (mass,sigma) values lay
  // ----------------------------------------------------
  // This must be done with respect to the width used in the computation of the probability distribution,
  // so that the bin 0 really matches the bin 0 of that distribution.
  //  double fracMass = (mass-(ResMass[iRes]-ResHalfWidth[iRes]))/(2*ResHalfWidth[iRes]);
  double fracMass = (mass - ResMinMass[iRes])/(2*ResHalfWidth[iRes]);
  if (debug>1) std::cout << std::setprecision(9)<<"mass ResMinMass[iRes] ResHalfWidth[iRes] ResHalfWidth[iRes]"
                    << mass << " "<<ResMinMass[iRes]<<" "<<ResHalfWidth[iRes]<<" "<<ResHalfWidth[iRes]<<std::endl;
  int iMassLeft  = (int)(fracMass*(double)nbins);
  int iMassRight = iMassLeft+1;
  double fracMassStep = (double)nbins*(fracMass - (double)iMassLeft/(double)nbins);
  if (debug>1) std::cout<<"nbins iMassLeft fracMass "<<nbins<<" "<<iMassLeft<<" "<<fracMass<<std::endl;

  // Simple protections for the time being: the region where we fit should not include
  // values outside the boundaries set by ResMass-ResHalfWidth : ResMass+ResHalfWidth
  // ---------------------------------------------------------------------------------
  if (iMassLeft<0) {
    edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassLeft="
                           << iMassLeft << "; mass = " << mass << " and bounds are " << ResMinMass[iRes]
                           << ":" << ResMinMass[iRes]+2*ResHalfWidth[iRes] << " - iMassLeft set to 0" << std::endl;
    iMassLeft  = 0;
    iMassRight = 1;
    insideProbMassWindow = false;
  }
  if (iMassRight>nbins) {
    edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassRight="
                           << iMassRight << "; mass = " << mass << " and bounds are " << ResMinMass[iRes]
                           << ":" << ResMass[iRes]+2*ResHalfWidth[iRes] << " - iMassRight set to " << nbins-1 << std::endl;
    iMassLeft  = nbins-1;
    iMassRight = nbins;
    insideProbMassWindow = false;
  }
  double fracSigma = (massResol/ResMaxSigma[iRes]);
  int iSigmaLeft = (int)(fracSigma*(double)nbins);
  int iSigmaRight = iSigmaLeft+1;
  double fracSigmaStep = (double)nbins * (fracSigma - (double)iSigmaLeft/(double)nbins);
  // std::cout << "massResol = " << massResol << std::endl;
  // std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
  // std::cout << "fracSigma = " << fracSigma << std::endl;
  // std::cout << "nbins = " << nbins << std::endl;
  // std::cout << "ISIGMALEFT = " << iSigmaLeft << std::endl;
  // std::cout << "ISIGMARIGHT = " << iSigmaRight << std::endl;
  // std::cout << "fracSigmaStep = " << fracSigmaStep << std::endl;

  // Simple protections for the time being: they should not affect convergence, since
  // ResMaxSigma is set to very large values, and if massResol exceeds them the fit
  // should not get any prize for that (for large sigma, the prob. distr. becomes flat)
  // ----------------------------------------------------------------------------------
  if (iSigmaLeft<0) {
    edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaLeft="
                           << iSigmaLeft << ", with massResol = " << massResol << " and ResMaxSigma[iRes] = "
                           << ResMaxSigma[iRes] << " -  iSigmaLeft set to 0" << std::endl;
    iSigmaLeft  = 0;
    iSigmaRight = 1;
  }
  if (iSigmaRight>nbins ) {
    if (counter_resprob<100)
      edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaRight="
                             << iSigmaRight << ", with massResol = " << massResol << " and ResMaxSigma[iRes] = "
                             << ResMaxSigma[iRes] << " -  iSigmaRight set to " << nbins-1 << std::endl;
    iSigmaLeft  = nbins-1;
    iSigmaRight = nbins;
  }

  // If f11,f12,f21,f22 are the values at the four corners, one finds by linear interpolation the
  // formula below for PS
  // --------------------------------------------------------------------------------------------
  if( insideProbMassWindow ) {
    double f11 = 0.;
    if (GLnorm[iY][iSigmaLeft]>0)
      f11 = GLvalue[iY][iMassLeft][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
    double f12 = 0.;
    if (GLnorm[iY][iSigmaRight]>0)
      f12 = GLvalue[iY][iMassLeft][iSigmaRight] / GLnorm[iY][iSigmaRight];
    double f21 = 0.;
    if (GLnorm[iY][iSigmaLeft]>0)
      f21 = GLvalue[iY][iMassRight][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
    double f22 = 0.;
    if (GLnorm[iY][iSigmaRight]>0)
      f22 = GLvalue[iY][iMassRight][iSigmaRight] / GLnorm[iY][iSigmaRight];
    PS = f11 + (f12-f11)*fracSigmaStep + (f21-f11)*fracMassStep +
      (f22-f21-f12+f11)*fracMassStep*fracSigmaStep;
    if (PS>0.1 || debug>1) LogDebug("MuScleFitUtils") << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
                                                      << f11 << " " << f12 << " " << f21 << " " << f22 << " "
                                                      << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
                                                      << iSigmaLeft << " " << iSigmaRight
                                                      << " GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
                                                      << " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;

//     if (PS>0.1) std::cout << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
//                      << f11 << " " << f12 << " " << f21 << " " << f22 << " "
//                      << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
//                      << iSigmaLeft << " " << iSigmaRight << " GLV,GLN="
//                      << GLvalue[iY][iMassLeft][iSigmaLeft]
//                      << " " << GLnorm[iY][iSigmaLeft] << std::endl;

  }
  else {
    edm::LogInfo("probability") << "outside mass probability window. Setting PS["<<iRes<<"] = 0" << std::endl;
  }

//   if( PS != PS ) {
//     std::cout << "ERROR: PS = " << PS << " for iRes = " << iRes << std::endl;

//     std::cout << "mass = " << mass << ", massResol = " << massResol << std::endl;
//     std::cout << "fracMass = " << fracMass << ", iMassLeft = " << iMassLeft
//          << ", iMassRight = " << iMassRight << ", fracMassStep = " << fracMassStep << std::endl;
//     std::cout << "fracSigma = " << fracSigma << ", iSigmaLeft = " << iSigmaLeft
//          << ", iSigmaRight = " << iSigmaRight << ", fracSigmaStep = " << fracSigmaStep << std::endl;
//     std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;

//     std::cout << "GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
//          << " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;
//   }

  return PS;
}