#include <CrossSectionHandler.h>

Public Member Functions
void	addParameters (std::vector< double > &initpar)
	Inputs the vars in a vector. More...

	CrossSectionHandler (const std::vector< double > &crossSection, const std::vector< int > &resfind)

unsigned int	parNum ()

std::vector< double >	relativeCrossSections (const double *variables, const std::vector< int > &resfind)
	Perform a variable transformation from N-1 to relative cross sections. More...

bool	releaseParameters (TMinuit &rmin, const std::vector< int > &resfind, const std::vector< int > &parfix, const int *ind, const int iorder, const unsigned int shift)
	Use the information in resfind, parorder and parfix to release the N-1 variables. More...

void	setParameters (double Start, double Step, double Mini, double Maxi, int ind, TString parname, const std::vector< double > &parCrossSection, const std::vector< int > &parCrossSectionOrder, const std::vector< int > &resfind)
	Initializes the arrays needed by Minuit. More...

Protected Member Functions
void	computeRelativeCrossSections (const std::vector< double > &crossSection, const std::vector< int > &resfind)

void	imposeConstraint ()
	Change of variables so that we move from N to N-1 variables using the constraint that Sum(x_i) = 1. More...

Protected Attributes
unsigned int	numberOfResonances_

unsigned int	parNum_

std::vector< double >	relativeCrossSectionVec_

std::vector< double >	vars_

Friends
class	TestCrossSectionHandler

Detailed Description

Definition at line 27 of file CrossSectionHandler.h.

Constructor & Destructor Documentation

CrossSectionHandler::CrossSectionHandler	(	const std::vector< double > &	crossSection,
		const std::vector< int > &	resfind
	)

inline

Definition at line 33 of file CrossSectionHandler.h.

References computeRelativeCrossSections(), imposeConstraint(), parNum_, and vars_.

                                                                                               :
     parNum_(0),
     numberOfResonances_(resfind.size())
   {
     // The number of parameters is the number of fitted resonances minus 1
     std::vector<int>::const_iterator it = resfind.begin();
     for( ; it != resfind.end(); ++it ) {
       if( *it != 0 ) ++parNum_;
     }
     if( parNum_ > 0 ) parNum_ = parNum_ - 1;
 
     vars_.resize(parNum_);
 
     computeRelativeCrossSections(crossSection, resfind);
     imposeConstraint();
   }

Member Function Documentation

void CrossSectionHandler::addParameters ( std::vector< double > & initpar )

inline

Inputs the vars in a vector.

Definition at line 51 of file CrossSectionHandler.h.

References vars_.

Referenced by MuScleFit::duringFastLoop().

   {
     std::vector<double>::const_iterator it = vars_.begin();
     for( ; it != vars_.end(); ++it ) {
       initpar.push_back(*it);
     }
   }

void CrossSectionHandler::computeRelativeCrossSections	(	const std::vector< double > &	crossSection,
		const std::vector< int > &	resfind
	)

inlineprotected

Initializes the relative cross sections for the range of resonances in [minRes, maxRes]. (note that both minRes and maxRes are included).
Also sets the lock on resonances. If only one of the resonances in the range is fitted its relative cross section will be 1 and it will not be fitted. If there are more than one only those that are fitted will have the relative cross section parameters unlocked during the fit.

Definition at line 166 of file CrossSectionHandler.h.

References ires, PostProcessor_cff::normalization, and relativeCrossSectionVec_.

Referenced by CrossSectionHandler(), and setParameters().

   {
     relativeCrossSectionVec_.clear();
     double normalization = 0.;
     for( unsigned int ires = 0; ires < resfind.size(); ++ires ) {
       if( resfind[ires] ) {
         normalization += crossSection[ires];
       }
     }
     if( normalization != 0. ) {
       for( unsigned int ires = 0; ires < resfind.size(); ++ires ) {
         if( resfind[ires] ) {
           relativeCrossSectionVec_.push_back(crossSection[ires]/normalization);
         }
       }
     }
   }

void CrossSectionHandler::imposeConstraint ( )

inlineprotected

Change of variables so that we move from N to N-1 variables using the constraint that Sum(x_i) = 1.

Definition at line 185 of file CrossSectionHandler.h.

References parNum_, relativeCrossSectionVec_, and vars_.

Referenced by CrossSectionHandler(), and setParameters().

   {
     if( parNum_ > 0 ) {
       for( unsigned int iVar = 0; iVar < parNum_; ++iVar ) {
         vars_[iVar] = relativeCrossSectionVec_[iVar+1]/relativeCrossSectionVec_[iVar];
       }
     }
   }

unsigned int CrossSectionHandler::parNum ( )

inline

Definition at line 115 of file CrossSectionHandler.h.

References parNum_.

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

   {
     return parNum_;
   }

std::vector<double> CrossSectionHandler::relativeCrossSections	(	const double *	variables,
		const std::vector< int > &	resfind
	)

inline

Perform a variable transformation from N-1 to relative cross sections.

Definition at line 121 of file CrossSectionHandler.h.

References mps_fire::i, numberOfResonances_, parNum_, and relativeCrossSectionVec_.

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

   {
     // parNum_ is 0 in two cases:
     // 1) if only one resonance is being fitted, in which case the relative cross section is
     // fixed to one and there is no need to recompute it
     // 2) if no resonance is being fitted, in which case all the relative cross sections will
     // be set to 0.
     // In both cases there is no need to make the transformation of variables.
     if( parNum_ != 0 ) {
       double * partialProduct = new double[numberOfResonances_];
       double norm = 0.;
       // Loop on all relative cross sections (that are parNum_+1)
       for( unsigned int i=0; i<parNum_+1; ++i ) {
         partialProduct[i] = std::accumulate(variables, variables + i, 1., std::multiplies<double>());
         norm += partialProduct[i];
       }
       for( unsigned int i=0; i<parNum_+1; ++i ) {
         relativeCrossSectionVec_[i] = partialProduct[i]/norm;
       }
       delete[] partialProduct;
     }
 
     std::vector<double> allRelativeCrossSections;
     std::vector<int>::const_iterator it = resfind.begin();
     int smallerVectorIndex = 0;
     for( ; it != resfind.end(); ++it ) {
       if( *it == 0 ) {
         allRelativeCrossSections.push_back(0.);
       }
       else {
         allRelativeCrossSections.push_back(relativeCrossSectionVec_[smallerVectorIndex]);
         ++smallerVectorIndex;
       }
     }
 
     return allRelativeCrossSections;
   }

bool CrossSectionHandler::releaseParameters	(	TMinuit &	rmin,
		const std::vector< int > &	resfind,
		const std::vector< int > &	parfix,
		const int *	ind,
		const int	iorder,
		const unsigned int	shift
	)

inline

Use the information in resfind, parorder and parfix to release the N-1 variables.

Definition at line 94 of file CrossSectionHandler.h.

References mps_fire::i, and numberOfResonances_.

Referenced by MuScleFitUtils::minimizeLikelihood().

   {
     // Find the number of free cross section parameters in this iteration
     unsigned int freeParNum = 0;
     for( unsigned int ipar=0; ipar<numberOfResonances_; ++ipar ) {
       if( (parfix[shift+ipar]==0) && (ind[shift+ipar]<=iorder) && (resfind[ipar] == 1) ) {
         ++freeParNum;
       }
     }
     if( freeParNum > 0 ) {
       freeParNum = freeParNum - 1;
       // Free only the first (freeParNum-1) of the N-1 variables
       for( unsigned int i=0; i<freeParNum; ++i ) {
         rmin.Release( shift+i );
       }
       return true;
     }
     return false;
   }

void CrossSectionHandler::setParameters	(	double *	Start,
		double *	Step,
		double *	Mini,
		double *	Maxi,
		int *	ind,
		TString *	parname,
		const std::vector< double > &	parCrossSection,
		const std::vector< int > &	parCrossSectionOrder,
		const std::vector< int > &	resfind
	)

inline

Initializes the arrays needed by Minuit.

Definition at line 60 of file CrossSectionHandler.h.

References computeRelativeCrossSections(), imposeConstraint(), numberOfResonances_, parNum_, and vars_.

Referenced by MuScleFitUtils::minimizeLikelihood().

   {
     computeRelativeCrossSections(parCrossSection, resfind);
     imposeConstraint();
 
     double thisStep[] = {0.001, 0.001, 0.001, 0.001, 0.001};
     TString thisParName[] = {"cross section var 1",
                              "cross section var 2",
                              "cross section var 3",
                              "cross section var 4",
                              "cross section var 5"};
     double thisMini[] = {0., 0., 0., 0., 0.};
     double thisMaxi[] = {1000., 1000., 1000., 1000., 1000.};
 
     // This is used to unlock the parameters in a given order. It is not
     // a TMinuit parameter, but a MuScleFit parameter.
     for( unsigned int iPar=0; iPar<numberOfResonances_; ++iPar ) {
       ind[iPar] = parCrossSectionOrder[iPar];
     }
 
     if( parNum_ > 0 ) {
       for( unsigned int iPar=0; iPar<parNum_; ++iPar ) {
         Start[iPar] = vars_[iPar];
         Step[iPar] = thisStep[iPar];
         Mini[iPar] = thisMini[iPar];
         Maxi[iPar] = thisMaxi[iPar];
         parname[iPar] = thisParName[iPar];
       }
     }
   }