MomentumScaleCalibration/interface/CrossSectionHandler.h

0001 /**
0002  * Handles the cross sections for MuScleFit. <br>
0003  * What counts in the fit is the ratio of the cross sections. However it depends on which resonances are used in the fit.
0004  * If we are fitting only the Upsilon(1S), for example, we do not need to consider the cross section ratio, because the probability
0005  * of the other resonances will be 0. This is useful when running on MC to test the algorithm. <br>
0006  *
0007  * The constructor receives the array of cross sections and the array of resfind that tells which of the resonances are considered in the fit. <br>
0008  * It builds the relative cross section parameters for each of the resonance and it has a method that unlocks the parameter accordingly. <br>
0009  * If for example only the Upsilon(1S) is fitted, the relative cross section will be 1 and it will remain fixed. <br>
0010  * The relative cross sections are fitted only when a background fit is done. <br>
0011  *
0012  * Note that this handles only the initialization of the cross sections, so that it is consistent with the fitted resonances, and
0013  * the fix/release of the cross section parameters. <br>
0014  *
0015  * This assumes that resfind is the same during all the processing (it is saved internally when received in the constructor).
0016  */
0017
0018 #include <vector>
0019 #include <numeric>
0020
0021 #include "TString.h"
0022 #include "TMinuit.h"
0023
0024 // Unit test for testing CrossSectionHandler
0025 class TestCrossSectionHandler;
0026
0027 class CrossSectionHandler {
0028   // For tests
0029   friend class TestCrossSectionHandler;
0030
0031 public:
0032   CrossSectionHandler(const std::vector<double>& crossSection, const std::vector<int>& resfind)
0033       : parNum_(0), numberOfResonances_(resfind.size()) {
0034     // The number of parameters is the number of fitted resonances minus 1
0035     std::vector<int>::const_iterator it = resfind.begin();
0036     for (; it != resfind.end(); ++it) {
0037       if (*it != 0)
0038         ++parNum_;
0039     }
0040     if (parNum_ > 0)
0041       parNum_ = parNum_ - 1;
0042
0043     vars_.resize(parNum_);
0044
0045     computeRelativeCrossSections(crossSection, resfind);
0046     imposeConstraint();
0047   }
0048
0049   /// Inputs the vars in a vector
0050   void addParameters(std::vector<double>& initpar) {
0051     std::vector<double>::const_iterator it = vars_.begin();
0052     for (; it != vars_.end(); ++it) {
0053       initpar.push_back(*it);
0054     }
0055   }
0056
0057   /// Initializes the arrays needed by Minuit
0058   void setParameters(double* Start,
0059                      double* Step,
0060                      double* Mini,
0061                      double* Maxi,
0062                      int* ind,
0063                      TString* parname,
0064                      const std::vector<double>& parCrossSection,
0065                      const std::vector<int>& parCrossSectionOrder,
0066                      const std::vector<int>& resfind) {
0067     computeRelativeCrossSections(parCrossSection, resfind);
0068     imposeConstraint();
0069
0070     double thisStep[] = {0.001, 0.001, 0.001, 0.001, 0.001};
0071     TString thisParName[] = {"cross section var 1",
0072                              "cross section var 2",
0073                              "cross section var 3",
0074                              "cross section var 4",
0075                              "cross section var 5"};
0076     double thisMini[] = {0., 0., 0., 0., 0.};
0077     double thisMaxi[] = {1000., 1000., 1000., 1000., 1000.};
0078
0079     // This is used to unlock the parameters in a given order. It is not
0080     // a TMinuit parameter, but a MuScleFit parameter.
0081     for (unsigned int iPar = 0; iPar < numberOfResonances_; ++iPar) {
0082       ind[iPar] = parCrossSectionOrder[iPar];
0083     }
0084
0085     if (parNum_ > 0) {
0086       for (unsigned int iPar = 0; iPar < parNum_; ++iPar) {
0087         Start[iPar] = vars_[iPar];
0088         Step[iPar] = thisStep[iPar];
0089         Mini[iPar] = thisMini[iPar];
0090         Maxi[iPar] = thisMaxi[iPar];
0091         parname[iPar] = thisParName[iPar];
0092       }
0093     }
0094   }
0095
0096   /// Use the information in resfind, parorder and parfix to release the N-1 variables
0097   bool releaseParameters(TMinuit& rmin,
0098                          const std::vector<int>& resfind,
0099                          const std::vector<int>& parfix,
0100                          const int* ind,
0101                          const int iorder,
0102                          const unsigned int shift) {
0103     // Find the number of free cross section parameters in this iteration
0104     unsigned int freeParNum = 0;
0105     for (unsigned int ipar = 0; ipar < numberOfResonances_; ++ipar) {
0106       if ((parfix[shift + ipar] == 0) && (ind[shift + ipar] <= iorder) && (resfind[ipar] == 1)) {
0107         ++freeParNum;
0108       }
0109     }
0110     if (freeParNum > 0) {
0111       freeParNum = freeParNum - 1;
0112       // Free only the first (freeParNum-1) of the N-1 variables
0113       for (unsigned int i = 0; i < freeParNum; ++i) {
0114         rmin.Release(shift + i);
0115       }
0116       return true;
0117     }
0118     return false;
0119   }
0120
0121   inline unsigned int parNum() { return parNum_; }
0122
0123   /// Perform a variable transformation from N-1 to relative cross sections
0124   std::vector<double> relativeCrossSections(const double* variables, const std::vector<int>& resfind) {
0125     // parNum_ is 0 in two cases:
0126     // 1) if only one resonance is being fitted, in which case the relative cross section is
0127     // fixed to one and there is no need to recompute it
0128     // 2) if no resonance is being fitted, in which case all the relative cross sections will
0129     // be set to 0.
0130     // In both cases there is no need to make the transformation of variables.
0131     if (parNum_ != 0) {
0132       double* partialProduct = new double[numberOfResonances_];
0133       double norm = 0.;
0134       // Loop on all relative cross sections (that are parNum_+1)
0135       for (unsigned int i = 0; i < parNum_ + 1; ++i) {
0136         partialProduct[i] = std::accumulate(variables, variables + i, 1., std::multiplies<double>());
0137         norm += partialProduct[i];
0138       }
0139       for (unsigned int i = 0; i < parNum_ + 1; ++i) {
0140         relativeCrossSectionVec_[i] = partialProduct[i] / norm;
0141       }
0142       delete[] partialProduct;
0143     }
0144
0145     std::vector<double> allRelativeCrossSections;
0146     std::vector<int>::const_iterator it = resfind.begin();
0147     int smallerVectorIndex = 0;
0148     for (; it != resfind.end(); ++it) {
0149       if (*it == 0) {
0150         allRelativeCrossSections.push_back(0.);
0151       } else {
0152         allRelativeCrossSections.push_back(relativeCrossSectionVec_[smallerVectorIndex]);
0153         ++smallerVectorIndex;
0154       }
0155     }
0156
0157     return allRelativeCrossSections;
0158   }
0159
0160 protected:
0161   /**
0162    * Initializes the relative cross sections for the range of resonances in [minRes, maxRes]. (note that both minRes and maxRes are included). <br>
0163    * 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
0164    * be fitted. If there are more than one only those that are fitted will have the relative cross section parameters unlocked during the fit.
0165    */
0166   // void computeRelativeCrossSections(const double crossSection[], const std::vector<int> resfind, const unsigned int minRes, const unsigned int maxRes)
0167   void computeRelativeCrossSections(const std::vector<double>& crossSection, const std::vector<int>& resfind) {
0168     relativeCrossSectionVec_.clear();
0169     double normalization = 0.;
0170     for (unsigned int ires = 0; ires < resfind.size(); ++ires) {
0171       if (resfind[ires]) {
0172         normalization += crossSection[ires];
0173       }
0174     }
0175     if (normalization != 0.) {
0176       for (unsigned int ires = 0; ires < resfind.size(); ++ires) {
0177         if (resfind[ires]) {
0178           relativeCrossSectionVec_.push_back(crossSection[ires] / normalization);
0179         }
0180       }
0181     }
0182   }
0183
0184   /// Change of variables so that we move from N to N-1 variables using the constraint that Sum(x_i) = 1.
0185   void imposeConstraint() {
0186     if (parNum_ > 0) {
0187       for (unsigned int iVar = 0; iVar < parNum_; ++iVar) {
0188         vars_[iVar] = relativeCrossSectionVec_[iVar + 1] / relativeCrossSectionVec_[iVar];
0189       }
0190     }
0191   }
0192
0193   // Data members
0194   std::vector<double> relativeCrossSectionVec_;
0195   std::vector<double> vars_;
0196   unsigned int parNum_;
0197   unsigned int numberOfResonances_;
0198 };