/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/MuonAnalysis/MomentumScaleCalibration/interface/CrossSectionHandler.h

Go to the documentation of this file.

#include <vector>
#include <numeric>

#include "TString.h"
#include "TMinuit.h"

class CrossSectionHandler
{
public:
  CrossSectionHandler(const std::vector<double> & crossSection, const std::vector<int> & resfind) :
    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();
  }

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

  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 )
  {
    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];
      }
    }
  }

  bool releaseParameters( TMinuit & rmin, const std::vector<int> & resfind, const std::vector<int> & parfix,
                          const int * ind, const int iorder, const unsigned int shift )
  {
    // 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;
  }

  inline unsigned int parNum()
  {
    return parNum_;
  }

  std::vector<double> relativeCrossSections( const double * variables, const std::vector<int> & resfind )
  {
    // 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;
  }

protected:
  // void computeRelativeCrossSections(const double crossSection[], const std::vector<int> resfind, const unsigned int minRes, const unsigned int maxRes)
  void computeRelativeCrossSections(const std::vector<double> & crossSection, const std::vector<int> & resfind)
  {
    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 imposeConstraint()
  {
    if( parNum_ > 0 ) {
      for( unsigned int iVar = 0; iVar < parNum_; ++iVar ) {
        vars_[iVar] = relativeCrossSectionVec_[iVar+1]/relativeCrossSectionVec_[iVar];
      }
    }
  }

  // Data members
  std::vector<double> relativeCrossSectionVec_;
  std::vector<double> vars_;
  unsigned int parNum_;
  unsigned int numberOfResonances_;
};