#include <BackgroundHandler.h>

Public Member Functions
std::pair< double, double >	backgroundFunction (const bool doBackgroundFit, const double parval, const int resTotNum, const int ires, const bool resConsidered, const double *ResMass, const double ResHalfWidth[], const int MuonType, const double &mass, const double &eta1, const double &eta2)

	BackgroundHandler (const std::vector< int > &identifiers, const std::vector< double > &leftWindowBorders, const std::vector< double > &rightWindowBorders, const double ResMass, const double massWindowHalfWidth)

bool	checkBackgroundWindow (const double &mass, const int iRegion)
	Check if the mass value is inside the given background region. More...

void	countEventsInAllWindows (const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &muonPairs, const double &weight)

void	initializeParNums ()
	Initialize the parNums to be used in the shifts of parval. More...

int	regionsParNum ()
	Returns the total number of parameters used for the regions. More...

void	rescale (std::vector< double > &parBgr, const double ResMass, const double massWindowHalfWidth, const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &muonPairs, const double &weight=1.)

double	resMass (const bool doBackgroundFit, const int ires)

void	setParameters (double Start, double Step, double Mini, double Maxi, int ind, TString parname, const std::vector< double > &parBgr, const std::vector< int > &parBgrOrder, const int muonType)
	Sets initial parameters for all the functions. More...

bool	unlockParameter (const std::vector< int > &resfind, const unsigned int ipar)
	returns true if the parameter is to be unlocked More...

std::pair< double, double >	windowBorders (const bool doBackgroundFit, const int ires)
	Returns the appropriate window borders depending on whether the background is being fitted and on the resonance. More...

	~BackgroundHandler ()

Private Member Functions
void	consistencyCheck (const std::vector< int > &identifiers, const std::vector< double > &leftWindowBorders, const std::vector< double > &rightWindowBorders) const throw (cms::Exception)
	Used to check the consistency of passed parameters. More...

Private Attributes
std::vector< MassWindow >	backgroundWindow_

int	parNumsRegions_ [3]

int	parNumsResonances_ [6]

int	regToResHW_ [3]

std::vector< MassWindow >	resonanceWindow_

int	resToReg_ [6]

Detailed Description

This class is used to handle the different background functions for the different regions.
It uses the backgroundFunctions defined in Functions.h and the backgroundFunctionService defined in Functions.cc.
More details are in the description of backgroundFunctionBase in Functions.h.

A bool selects if to use the regions functions or the resonances functions.

Definition at line 29 of file BackgroundHandler.h.

Constructor & Destructor Documentation

BackgroundHandler::BackgroundHandler	(	const std::vector< int > &	identifiers,
		const std::vector< double > &	leftWindowBorders,
		const std::vector< double > &	rightWindowBorders,
		const double *	ResMass,
		const double *	massWindowHalfWidth
	)

Definition at line 12 of file BackgroundHandler.cc.

References backgroundFunctionService(), backgroundWindow_, consistencyCheck(), i, cmsHarvester::index, initializeParNums(), regToResHW_, resonanceWindow_, and resToReg_.

 {
   // : leftWindowFactors_(leftWindowFactors),
   // rightWindowFactors_(rightWindowFactors)
 
   // Define the correspondence between regions and halfWidth to use
   // Defines also the function type to use (but they are checked to be consistent over a region)
   regToResHW_[0] = 0; // Region 0 use the one from Z
   regToResHW_[1] = 3; // Region 1 use the one from Upsilon1S
   regToResHW_[2] = 5; // Region 2 use the one from J/Psi
 
   // Define the correspondence between resonances and regions
   resToReg_[0] = 0; // Z
   resToReg_[1] = 1; // Upsilon3S
   resToReg_[2] = 1; // Upsilon2S
   resToReg_[3] = 1; // Upsilon1S
   resToReg_[4] = 2; // Psi2S
   resToReg_[5] = 2; // J/Psi
 
   // Throws cms::Exception("Configuration") in case the parameters are not what is expected
   consistencyCheck(identifiers, leftWindowBorders, rightWindowBorders);
 
   // Build the resonance windows
   for( unsigned int i=0; i<6; ++i ) {
     double mass = ResMass[i];
     double lowerLimit = mass - massWindowHalfWidth[i];
     double upperLimit = mass + massWindowHalfWidth[i];
     resonanceWindow_.push_back( MassWindow( mass, lowerLimit, upperLimit,
                                             std::vector<unsigned int>(1,i),
                                             backgroundFunctionService(identifiers[resToReg_[i]],
                                                                       lowerLimit,
                                                                       upperLimit) ) );
   }
 
   // Build the background windows
   // ----------------------------
   // Compute the mass center of each region
   double resMassForRegion[3];
   resMassForRegion[0] = ResMass[0];
   resMassForRegion[1] = (ResMass[1]+ResMass[2]+ResMass[3])/3;
   resMassForRegion[2] = (ResMass[4]+ResMass[5])/2;
 
   // Define which resonance is in which background window
   std::vector<std::vector<unsigned int> > indexes;
   indexes.push_back(std::vector<unsigned int>(1,0));
   indexes.push_back(std::vector<unsigned int>());
   for( int i=1; i<=3; ++i ) { indexes[1].push_back(i); }
   indexes.push_back(std::vector<unsigned int>());
   for( int i=4; i<=5; ++i ) { indexes[2].push_back(i); }
 
   unsigned int i=0;
   typedef std::vector<unsigned int> indexVec;
   BOOST_FOREACH(const indexVec & index, indexes) {
     //     double lowerLimit = resMassForRegion[i] - massWindowHalfWidth[regToResHW_[i]]*leftWindowFactors[i];
     //     double upperLimit = resMassForRegion[i] + massWindowHalfWidth[regToResHW_[i]]*rightWindowFactors[i];
     //     backgroundWindow_.push_back( MassWindow( resMassForRegion[i], lowerLimit, upperLimit, index,
     //                                              backgroundFunctionService(identifiers[i], lowerLimit, upperLimit ) ) );
     backgroundWindow_.push_back( MassWindow( resMassForRegion[i], leftWindowBorders[i], rightWindowBorders[i], index,
                                              backgroundFunctionService(identifiers[i], leftWindowBorders[i], rightWindowBorders[i] ) ) );
     ++i;
   }
   // Initialize the parNums to be used in the shifts of parval
   initializeParNums();
 }

BackgroundHandler::~BackgroundHandler ( )

Definition at line 81 of file BackgroundHandler.cc.

82 {

83 }

Member Function Documentation

std::pair< double, double > BackgroundHandler::backgroundFunction	(	const bool	doBackgroundFit,
		const double *	parval,
		const int	resTotNum,
		const int	ires,
		const bool *	resConsidered,
		const double *	ResMass,
		const double	ResHalfWidth[],
		const int	MuonType,
		const double &	mass,
		const double &	eta1,
		const double &	eta2
	)

Returns the background fraction parameter (parBgr[0], but shifted to the correct function) and the value returned by the background function.
Depending on the value of doBackgroundFit it returns the values for the regions or the resonances.

Definition at line 226 of file BackgroundHandler.cc.

References backgroundWindow_, ires, parNumsRegions_, parNumsResonances_, resonanceWindow_, and resToReg_.

Referenced by MuScleFitUtils::massProb().

 {
   if( doBackgroundFit ) {
     // Return the values for the region
     int iReg = resToReg_[ires];
     // return std::make_pair( parval[parNumsRegions_[iReg]] * backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), resEta),
     return std::make_pair( parval[parNumsRegions_[iReg]] * backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), eta1, eta2 ),
                    (*(backgroundWindow_[iReg].backgroundFunction()))( &(parval[parNumsRegions_[iReg]]), mass, eta1, eta2 ) );
     // return std::make_pair( backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), eta1, eta2),
     //             (*(backgroundWindow_[iReg].backgroundFunction()))( &(parval[parNumsRegions_[iReg]]), mass, eta1, eta2 ) );
   }
   // Return the values for the resonance
   // return std::make_pair( parval[parNumsResonances_[ires]] * resonanceWindow_[ires].backgroundFunction()->fracVsEta(&(parval[parNumsResonances_[ires]]), resEta),
   //             (*(resonanceWindow_[ires].backgroundFunction()))( &(parval[parNumsResonances_[ires]]), mass, resEta ) );
   return std::make_pair( parval[parNumsResonances_[ires]] * resonanceWindow_[ires].backgroundFunction()->fracVsEta(&(parval[parNumsResonances_[ires]]), eta1, eta2),
              (*(resonanceWindow_[ires].backgroundFunction()))( &(parval[parNumsResonances_[ires]]), mass, eta1, eta2 ) );
 }

bool BackgroundHandler::checkBackgroundWindow	(	const double &	mass,
		const int	iRegion
	)

inline

Check if the mass value is inside the given background region.

Definition at line 49 of file BackgroundHandler.h.

References backgroundWindow_.

   {
     return backgroundWindow_[iRegion].isIn(mass);
   }

void BackgroundHandler::consistencyCheck	(	const std::vector< int > &	identifiers,
		const std::vector< double > &	leftWindowBorders,
		const std::vector< double > &	rightWindowBorders
	)		const
throw	(	cms::Exception
	)

private

Used to check the consistency of passed parameters.

Definition at line 274 of file BackgroundHandler.cc.

References Exception.

Referenced by BackgroundHandler().

 {
   if( leftWindowBorders.size() != rightWindowBorders.size() ) {
     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = " << leftWindowBorders.size()
                                           << " != rightWindowBorders.size() = " << rightWindowBorders.size() << std::endl;
   }
   if( leftWindowBorders.size() != 3 ) {
     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = rightWindowBorders.size() = "
                                           << leftWindowBorders.size() << " != 3" << std::endl;
   }
   if( identifiers.size() != 3 ) {
     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: identifiers must match the number of regions = 3" << std::endl;
   }
 }

void BackgroundHandler::countEventsInAllWindows	(	const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &	muonPairs,
		const double &	weight
	)

Definition at line 248 of file BackgroundHandler.cc.

References backgroundWindow_, MassWindow::count(), MassWindow::resetCounter(), resonanceWindow_, and histoStyle::weight.

Referenced by rescale().

 {
   // First reset all the counters
   BOOST_FOREACH(MassWindow & resonanceWindow, resonanceWindow_) {
     resonanceWindow.resetCounter();
   }
   // Count events in background windows
   BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_) {
     backgroundWindow.resetCounter();
   }
 
   // Now count the events in each window
   std::pair<lorentzVector,lorentzVector> muonPair;
   BOOST_FOREACH(muonPair, muonPairs) {
     // Count events in resonance windows
     BOOST_FOREACH(MassWindow & resonanceWindow, resonanceWindow_) {
       resonanceWindow.count((muonPair.first + muonPair.second).mass(), weight);
     }
     // Count events in background windows
     BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_) {
       backgroundWindow.count((muonPair.first + muonPair.second).mass(), weight);
     }
   }
 }

void BackgroundHandler::initializeParNums ( )

Initialize the parNums to be used in the shifts of parval.

Definition at line 85 of file BackgroundHandler.cc.

References backgroundWindow_, i, parNumsRegions_, parNumsResonances_, and resonanceWindow_.

Referenced by BackgroundHandler().

 {
   // Initialize the parNums to be used in the shifts of parval
   parNumsRegions_[0] = 0;
   for( unsigned int i=1; i<backgroundWindow_.size() ; ++i ) {
     parNumsRegions_[i] = parNumsRegions_[i-1] + backgroundWindow_[i-1].backgroundFunction()->parNum();
   }
   parNumsResonances_[0] = parNumsRegions_[2]+backgroundWindow_[2].backgroundFunction()->parNum();
   for( unsigned int i=1; i<resonanceWindow_.size() ; ++i ) {
     parNumsResonances_[i] = parNumsResonances_[i-1] + resonanceWindow_[i-1].backgroundFunction()->parNum();
   }
 }

int BackgroundHandler::regionsParNum ( )

inline

Returns the total number of parameters used for the regions.

Definition at line 43 of file BackgroundHandler.h.

References parNumsResonances_.

Referenced by MuScleFitUtils::minimizeLikelihood().

   {
     return parNumsResonances_[0];
   }

void BackgroundHandler::rescale	(	std::vector< double > &	parBgr,
		const double *	ResMass,
		const double *	massWindowHalfWidth,
		const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &	muonPairs,
		const double &	weight = `1.`
	)

Computes the rescaled parameters from the regions functions to the resonances functions. It takes into account the difference in intervals and rescales the parameters so that the fraction of events is correctly accounted for.
It uses the list of all muon pairs to compute the number of events in each resonance window.

Definition at line 171 of file BackgroundHandler.cc.

References MassWindow::backgroundFunction(), backgroundWindow_, countEventsInAllWindows(), gather_cfg::cout, MassWindow::events(), backgroundFunctionBase::functionForIntegral(), cmsHarvester::index, MassWindow::indexes(), MassWindow::lowerBound(), parNumsRegions_, parNumsResonances_, resonanceWindow_, resToReg_, and MassWindow::upperBound().

Referenced by MuScleFitUtils::minimizeLikelihood().

 {
   countEventsInAllWindows(muonPairs, weight);
 
   // Loop on all regions and on all the resonances of each region and compute the background fraction
   // for each resonance window.
   unsigned int iRegion = 0;
   BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_)
   {
     // Iterator pointing to the first parameter of this background function in the full set of parameters
     std::vector<double>::const_iterator parBgrIt = (parBgr.begin()+parNumsRegions_[iRegion]);
     TF1 * backgroundFunctionForIntegral = backgroundWindow.backgroundFunction()->functionForIntegral(parBgrIt);
     // WARNING: this expects the background fraction parameter to be parBgr[0] for all the background functions.
     double kOld = *parBgrIt;
     double Nbw = backgroundWindow.events();
     double Ibw = backgroundFunctionForIntegral->Integral(backgroundWindow.lowerBound(),
                                                          backgroundWindow.upperBound());
 
     // index is the index of the resonance in the background window
     BOOST_FOREACH( unsigned int index, *(backgroundWindow.indexes()) )
     {
       // First set all parameters of the resonance window background function to those of the corresponding region
       for( int iPar = 0; iPar < resonanceWindow_[index].backgroundFunction()->parNum(); ++iPar ) {
         parBgr[parNumsResonances_[index]+iPar] = parBgr[parNumsRegions_[resToReg_[index]]+iPar];
       }
       // Estimated fraction of events in the resonance window
       double Irw = backgroundFunctionForIntegral->Integral(resonanceWindow_[index].lowerBound(),
                                                            resonanceWindow_[index].upperBound());
       double Nrw = resonanceWindow_[index].events();
 
       // Ibw is 1 (to avoid effects from approximation errors we set it to 1 and do not include it in the computation).
       if( Nrw != 0 ) parBgr[parNumsResonances_[index]] = kOld*Nbw/Nrw*Irw;
       else parBgr[parNumsResonances_[index]] = 0.;
 
       // Protect against fluctuations of number of events which could cause the fraction to go above 1.
       if( parBgr[parNumsResonances_[index]] > 1. ) parBgr[parNumsResonances_[index]] = 1.;
 
       double kNew = parBgr[parNumsResonances_[index]];
       std::cout << "For resonance = " << index << std::endl;
       std::cout << "backgroundWindow.lowerBound = " << backgroundWindow.lowerBound() << std::endl;
       std::cout << "backgroundWindow.upperBound = " << backgroundWindow.upperBound() << std::endl;
       std::cout << "parNumsResonances_["<<index<<"] = " << parNumsResonances_[index] << std::endl;
       std::cout << "Nbw = " << Nbw << ", Ibw = " << Ibw << std::endl;
       std::cout << "Nrw = " << Nrw << ", Irw = " << Irw << std::endl;
       std::cout << "k = " << kOld << ", k' = " << parBgr[parNumsResonances_[index]] << std::endl;
       std::cout << "background fraction in background window = Nbw*k = " << Nbw*kOld << std::endl;
       std::cout << "background fraction in resonance window = Nrw*k' = " << Nrw*kNew << std::endl;
     }
     ++iRegion;
     delete backgroundFunctionForIntegral;
   }
 }

double BackgroundHandler::resMass	(	const bool	doBackgroundFit,
		const int	ires
	)

Returns the appropriate resMass value depending on whether the background is being fitted and on the resonance.
The resMass used for the region is the mean of the mass of the corresponding resonances, so for the Z is the same Z mass, for the Upsilons is the arithmetic mean of the Upsilon masses and the same for the J/Psi and Psi2S region.

Definition at line 159 of file BackgroundHandler.cc.

References backgroundWindow_, ires, resonanceWindow_, and resToReg_.

Referenced by MuScleFitUtils::massProb().

 {
   if( doBackgroundFit ) {
     // Fitting the background: use the regions
     return backgroundWindow_[resToReg_[ires]].mass();
   }
   else {
     // Not fitting the background: use the resonances
     return resonanceWindow_[ires].mass();
   }
 }

void BackgroundHandler::setParameters	(	double *	Start,
		double *	Step,
		double *	Mini,
		double *	Maxi,
		int *	ind,
		TString *	parname,
		const std::vector< double > &	parBgr,
		const std::vector< int > &	parBgrOrder,
		const int	muonType
	)

Sets initial parameters for all the functions.

Definition at line 98 of file BackgroundHandler.cc.

References backgroundWindow_, parNumsRegions_, parNumsResonances_, resonanceWindow_, and edm::shift.

Referenced by MuScleFitUtils::minimizeLikelihood().

 {
   std::vector<double>::const_iterator parBgrIt = parBgr.begin();
   std::vector<int>::const_iterator parBgrOrderIt = parBgrOrder.begin();
   // Set the parameters for the regions only if this is not a rescaling
   for( int iReg = 0; iReg < 3; ++iReg ) {
     int shift = parNumsRegions_[iReg];
     backgroundWindow_[iReg].backgroundFunction()->setParameters( &(Start[shift]), &(Step[shift]), &(Mini[shift]),
         &(Maxi[shift]), &(ind[shift]), &(parname[shift]),
         parBgrIt+shift, parBgrOrderIt+shift, muonType );
   }
   for( int iRes = 0; iRes < 6; ++iRes ) {
     // parNumsResonances is already shifted for the regions parameters
     int shift = parNumsResonances_[iRes];
     resonanceWindow_[iRes].backgroundFunction()->setParameters( &(Start[shift]), &(Step[shift]), &(Mini[shift]),
         &(Maxi[shift]), &(ind[shift]), &(parname[shift]),
         parBgrIt+shift, parBgrOrderIt+shift, muonType );
   }
 }

bool BackgroundHandler::unlockParameter	(	const std::vector< int > &	resfind,
		const unsigned int	ipar
	)

returns true if the parameter is to be unlocked

Definition at line 118 of file BackgroundHandler.cc.

References parNumsRegions_, and parNumsResonances_.

Referenced by MuScleFitUtils::minimizeLikelihood().

 {
   // parNumsRegions_ are shifted: [1] contains the number of parameters for 0 and so on.
   if( ipar < unsigned(parNumsRegions_[1]) && resfind[0] > 0 ) {
     return true;
   }
   if( ipar >= unsigned(parNumsRegions_[1]) && ipar < unsigned(parNumsRegions_[2]) && ( resfind[1] > 0 || resfind[2] > 0 || resfind[3] > 0 ) ) {
     return true;
   }
   // The first of parNumsResonances_ has the sum of parNums of the regions.
   if( ipar >= unsigned(parNumsRegions_[2]) && ipar < unsigned(parNumsResonances_[0]) && ( resfind[4] > 0 || resfind[5] > 0 ) ) {
     return true;
   }
   return false;
 }

std::pair< double, double > BackgroundHandler::windowBorders	(	const bool	doBackgroundFit,
		const int	ires
	)

Returns the appropriate window borders depending on whether the background is being fitted and on the resonance.

Definition at line 146 of file BackgroundHandler.cc.

References backgroundWindow_, resonanceWindow_, and resToReg_.

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

 {
   if( doBackgroundFit ) {
     // Fitting the background: use the regions
     return std::make_pair(backgroundWindow_[resToReg_[ires]].lowerBound(), backgroundWindow_[resToReg_[ires]].upperBound());
   }
   else {
     // Not fitting the background: use the resonances
     // return std::make_pair(1.,1.);
     return std::make_pair(resonanceWindow_[ires].lowerBound(), resonanceWindow_[ires].upperBound());
   }
 }