MomentumScaleCalibration/src/BackgroundHandler.cc

0001 #ifndef BackgroundHandler_cc
0002 #define BackgroundHandler_cc
0003
0004 #include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundHandler.h"
0005 #include <algorithm>
0006 #include <TF1.h>
0007 #include <iostream>
0008
0009 typedef reco::Particle::LorentzVector lorentzVector;
0010
0011 BackgroundHandler::BackgroundHandler(const std::vector<int>& identifiers,
0012                                      const std::vector<double>& leftWindowBorders,
0013                                      const std::vector<double>& rightWindowBorders,
0014                                      const double* ResMass,
0015                                      const double* massWindowHalfWidth) {
0016   // : leftWindowFactors_(leftWindowFactors),
0017   // rightWindowFactors_(rightWindowFactors)
0018
0019   // Define the correspondence between regions and halfWidth to use
0020   // Defines also the function type to use (but they are checked to be consistent over a region)
0021   regToResHW_[0] = 0;  // Region 0 use the one from Z
0022   regToResHW_[1] = 3;  // Region 1 use the one from Upsilon1S
0023   regToResHW_[2] = 5;  // Region 2 use the one from J/Psi
0024
0025   // Define the correspondence between resonances and regions
0026   resToReg_[0] = 0;  // Z
0027   resToReg_[1] = 1;  // Upsilon3S
0028   resToReg_[2] = 1;  // Upsilon2S
0029   resToReg_[3] = 1;  // Upsilon1S
0030   resToReg_[4] = 2;  // Psi2S
0031   resToReg_[5] = 2;  // J/Psi
0032
0033   // Throws cms::Exception("Configuration") in case the parameters are not what is expected
0034   consistencyCheck(identifiers, leftWindowBorders, rightWindowBorders);
0035
0036   // Build the resonance windows
0037   for (unsigned int i = 0; i < 6; ++i) {
0038     double mass = ResMass[i];
0039     double lowerLimit = mass - massWindowHalfWidth[i];
0040     double upperLimit = mass + massWindowHalfWidth[i];
0041     resonanceWindow_.push_back(
0042         MassWindow(mass,
0043                    lowerLimit,
0044                    upperLimit,
0045                    std::vector<unsigned int>(1, i),
0046                    backgroundFunctionService(identifiers[resToReg_[i]], lowerLimit, upperLimit)));
0047   }
0048
0049   // Build the background windows
0050   // ----------------------------
0051   // Compute the mass center of each region
0052   double resMassForRegion[3];
0053   resMassForRegion[0] = ResMass[0];
0054   resMassForRegion[1] = (ResMass[1] + ResMass[2] + ResMass[3]) / 3;
0055   resMassForRegion[2] = (ResMass[4] + ResMass[5]) / 2;
0056
0057   // Define which resonance is in which background window
0058   std::vector<std::vector<unsigned int> > indexes;
0059   indexes.push_back(std::vector<unsigned int>(1, 0));
0060   indexes.push_back(std::vector<unsigned int>());
0061   for (int i = 1; i <= 3; ++i) {
0062     indexes[1].push_back(i);
0063   }
0064   indexes.push_back(std::vector<unsigned int>());
0065   for (int i = 4; i <= 5; ++i) {
0066     indexes[2].push_back(i);
0067   }
0068
0069   unsigned int i = 0;
0070   typedef std::vector<unsigned int> indexVec;
0071   for (auto const& index : indexes) {
0072     //     double lowerLimit = resMassForRegion[i] - massWindowHalfWidth[regToResHW_[i]]*leftWindowFactors[i];
0073     //     double upperLimit = resMassForRegion[i] + massWindowHalfWidth[regToResHW_[i]]*rightWindowFactors[i];
0074     //     backgroundWindow_.push_back( MassWindow( resMassForRegion[i], lowerLimit, upperLimit, index,
0075     //                                              backgroundFunctionService(identifiers[i], lowerLimit, upperLimit ) ) );
0076     backgroundWindow_.push_back(
0077         MassWindow(resMassForRegion[i],
0078                    leftWindowBorders[i],
0079                    rightWindowBorders[i],
0080                    index,
0081                    backgroundFunctionService(identifiers[i], leftWindowBorders[i], rightWindowBorders[i])));
0082     ++i;
0083   }
0084   // Initialize the parNums to be used in the shifts of parval
0085   initializeParNums();
0086 }
0087
0088 BackgroundHandler::~BackgroundHandler() {}
0089
0090 void BackgroundHandler::initializeParNums() {
0091   // Initialize the parNums to be used in the shifts of parval
0092   parNumsRegions_[0] = 0;
0093   for (unsigned int i = 1; i < backgroundWindow_.size(); ++i) {
0094     parNumsRegions_[i] = parNumsRegions_[i - 1] + backgroundWindow_[i - 1].backgroundFunction()->parNum();
0095   }
0096   parNumsResonances_[0] = parNumsRegions_[2] + backgroundWindow_[2].backgroundFunction()->parNum();
0097   for (unsigned int i = 1; i < resonanceWindow_.size(); ++i) {
0098     parNumsResonances_[i] = parNumsResonances_[i - 1] + resonanceWindow_[i - 1].backgroundFunction()->parNum();
0099   }
0100 }
0101
0102 void BackgroundHandler::setParameters(double* Start,
0103                                       double* Step,
0104                                       double* Mini,
0105                                       double* Maxi,
0106                                       int* ind,
0107                                       TString* parname,
0108                                       const std::vector<double>& parBgr,
0109                                       const std::vector<int>& parBgrOrder,
0110                                       const int muonType) {
0111   std::vector<double>::const_iterator parBgrIt = parBgr.begin();
0112   std::vector<int>::const_iterator parBgrOrderIt = parBgrOrder.begin();
0113   // Set the parameters for the regions only if this is not a rescaling
0114   for (int iReg = 0; iReg < 3; ++iReg) {
0115     int shift = parNumsRegions_[iReg];
0116     backgroundWindow_[iReg].backgroundFunction()->setParameters(&(Start[shift]),
0117                                                                 &(Step[shift]),
0118                                                                 &(Mini[shift]),
0119                                                                 &(Maxi[shift]),
0120                                                                 &(ind[shift]),
0121                                                                 &(parname[shift]),
0122                                                                 parBgrIt + shift,
0123                                                                 parBgrOrderIt + shift,
0124                                                                 muonType);
0125   }
0126   for (int iRes = 0; iRes < 6; ++iRes) {
0127     // parNumsResonances is already shifted for the regions parameters
0128     int shift = parNumsResonances_[iRes];
0129     resonanceWindow_[iRes].backgroundFunction()->setParameters(&(Start[shift]),
0130                                                                &(Step[shift]),
0131                                                                &(Mini[shift]),
0132                                                                &(Maxi[shift]),
0133                                                                &(ind[shift]),
0134                                                                &(parname[shift]),
0135                                                                parBgrIt + shift,
0136                                                                parBgrOrderIt + shift,
0137                                                                muonType);
0138   }
0139 }
0140
0141 bool BackgroundHandler::unlockParameter(const std::vector<int>& resfind, const unsigned int ipar) {
0142   // parNumsRegions_ are shifted: [1] contains the number of parameters for 0 and so on.
0143   if (ipar < unsigned(parNumsRegions_[1]) && resfind[0] > 0) {
0144     return true;
0145   }
0146   if (ipar >= unsigned(parNumsRegions_[1]) && ipar < unsigned(parNumsRegions_[2]) &&
0147       (resfind[1] > 0 || resfind[2] > 0 || resfind[3] > 0)) {
0148     return true;
0149   }
0150   // The first of parNumsResonances_ has the sum of parNums of the regions.
0151   if (ipar >= unsigned(parNumsRegions_[2]) && ipar < unsigned(parNumsResonances_[0]) &&
0152       (resfind[4] > 0 || resfind[5] > 0)) {
0153     return true;
0154   }
0155   return false;
0156 }
0157
0158 // std::pair<double, double> BackgroundHandler::windowFactors( const bool doBackgroundFit, const int ires )
0159 // {
0160 //   if( doBackgroundFit ) {
0161 //     // Fitting the background: use the regions
0162 //     return std::make_pair(leftWindowFactors_[resToReg_[ires]], rightWindowFactors_[resToReg_[ires]]);
0163 //   }
0164 //   else {
0165 //     // Not fitting the background: use the resonances
0166 //     return std::make_pair(1.,1.);
0167 //   }
0168 // }
0169
0170 std::pair<double, double> BackgroundHandler::windowBorders(const bool doBackgroundFit, const int ires) {
0171   if (doBackgroundFit) {
0172     // Fitting the background: use the regions
0173     return std::make_pair(backgroundWindow_[resToReg_[ires]].lowerBound(),
0174                           backgroundWindow_[resToReg_[ires]].upperBound());
0175   } else {
0176     // Not fitting the background: use the resonances
0177     // return std::make_pair(1.,1.);
0178     return std::make_pair(resonanceWindow_[ires].lowerBound(), resonanceWindow_[ires].upperBound());
0179   }
0180 }
0181
0182 double BackgroundHandler::resMass(const bool doBackgroundFit, const int ires) {
0183   if (doBackgroundFit) {
0184     // Fitting the background: use the regions
0185     return backgroundWindow_[resToReg_[ires]].mass();
0186   } else {
0187     // Not fitting the background: use the resonances
0188     return resonanceWindow_[ires].mass();
0189   }
0190 }
0191
0192 void BackgroundHandler::rescale(
0193     std::vector<double>& parBgr,
0194     const double* ResMass,
0195     const double* massWindowHalfWidth,
0196     const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >& muonPairs,
0197     const double& weight) {
0198   countEventsInAllWindows(muonPairs, weight);
0199
0200   // Loop on all regions and on all the resonances of each region and compute the background fraction
0201   // for each resonance window.
0202   unsigned int iRegion = 0;
0203   for (auto const& backgroundWindow : backgroundWindow_) {
0204     // Iterator pointing to the first parameter of this background function in the full set of parameters
0205     std::vector<double>::const_iterator parBgrIt = (parBgr.begin() + parNumsRegions_[iRegion]);
0206     TF1* backgroundFunctionForIntegral = backgroundWindow.backgroundFunction()->functionForIntegral(parBgrIt);
0207     // WARNING: this expects the background fraction parameter to be parBgr[0] for all the background functions.
0208     double kOld = *parBgrIt;
0209     double Nbw = backgroundWindow.events();
0210     double Ibw = backgroundFunctionForIntegral->Integral(backgroundWindow.lowerBound(), backgroundWindow.upperBound());
0211
0212     // index is the index of the resonance in the background window
0213     for (unsigned int index : *(backgroundWindow.indexes())) {
0214       // First set all parameters of the resonance window background function to those of the corresponding region
0215       for (int iPar = 0; iPar < resonanceWindow_[index].backgroundFunction()->parNum(); ++iPar) {
0216         parBgr[parNumsResonances_[index] + iPar] = parBgr[parNumsRegions_[resToReg_[index]] + iPar];
0217       }
0218       // Estimated fraction of events in the resonance window
0219       double Irw = backgroundFunctionForIntegral->Integral(resonanceWindow_[index].lowerBound(),
0220                                                            resonanceWindow_[index].upperBound());
0221       double Nrw = resonanceWindow_[index].events();
0222
0223       // Ibw is 1 (to avoid effects from approximation errors we set it to 1 and do not include it in the computation).
0224       if (Nrw != 0)
0225         parBgr[parNumsResonances_[index]] = kOld * Nbw / Nrw * Irw;
0226       else
0227         parBgr[parNumsResonances_[index]] = 0.;
0228
0229       // Protect against fluctuations of number of events which could cause the fraction to go above 1.
0230       if (parBgr[parNumsResonances_[index]] > 1.)
0231         parBgr[parNumsResonances_[index]] = 1.;
0232
0233       double kNew = parBgr[parNumsResonances_[index]];
0234       std::cout << "For resonance = " << index << std::endl;
0235       std::cout << "backgroundWindow.lowerBound = " << backgroundWindow.lowerBound() << std::endl;
0236       std::cout << "backgroundWindow.upperBound = " << backgroundWindow.upperBound() << std::endl;
0237       std::cout << "parNumsResonances_[" << index << "] = " << parNumsResonances_[index] << std::endl;
0238       std::cout << "Nbw = " << Nbw << ", Ibw = " << Ibw << std::endl;
0239       std::cout << "Nrw = " << Nrw << ", Irw = " << Irw << std::endl;
0240       std::cout << "k = " << kOld << ", k' = " << parBgr[parNumsResonances_[index]] << std::endl;
0241       std::cout << "background fraction in background window = Nbw*k = " << Nbw * kOld << std::endl;
0242       std::cout << "background fraction in resonance window = Nrw*k' = " << Nrw * kNew << std::endl;
0243     }
0244     ++iRegion;
0245     delete backgroundFunctionForIntegral;
0246   }
0247 }
0248
0249 std::pair<double, double> BackgroundHandler::backgroundFunction(const bool doBackgroundFit,
0250                                                                 const double* parval,
0251                                                                 const int resTotNum,
0252                                                                 const int ires,
0253                                                                 const bool* resConsidered,
0254                                                                 const double* ResMass,
0255                                                                 const double ResHalfWidth[],
0256                                                                 const int MuonType,
0257                                                                 const double& mass,
0258                                                                 const double& eta1,
0259                                                                 const double& eta2) {
0260   if (doBackgroundFit) {
0261     // Return the values for the region
0262     int iReg = resToReg_[ires];
0263     // return std::make_pair( parval[parNumsRegions_[iReg]] * backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), resEta),
0264     return std::make_pair(
0265         parval[parNumsRegions_[iReg]] *
0266             backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), eta1, eta2),
0267         (*(backgroundWindow_[iReg].backgroundFunction()))(&(parval[parNumsRegions_[iReg]]), mass, eta1, eta2));
0268     // return std::make_pair( backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), eta1, eta2),
0269     //             (*(backgroundWindow_[iReg].backgroundFunction()))( &(parval[parNumsRegions_[iReg]]), mass, eta1, eta2 ) );
0270   }
0271   // Return the values for the resonance
0272   // return std::make_pair( parval[parNumsResonances_[ires]] * resonanceWindow_[ires].backgroundFunction()->fracVsEta(&(parval[parNumsResonances_[ires]]), resEta),
0273   //             (*(resonanceWindow_[ires].backgroundFunction()))( &(parval[parNumsResonances_[ires]]), mass, resEta ) );
0274   return std::make_pair(
0275       parval[parNumsResonances_[ires]] *
0276           resonanceWindow_[ires].backgroundFunction()->fracVsEta(&(parval[parNumsResonances_[ires]]), eta1, eta2),
0277       (*(resonanceWindow_[ires].backgroundFunction()))(&(parval[parNumsResonances_[ires]]), mass, eta1, eta2));
0278 }
0279
0280 void BackgroundHandler::countEventsInAllWindows(
0281     const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >& muonPairs,
0282     const double& weight) {
0283   // First reset all the counters
0284   for (auto& resonanceWindow : resonanceWindow_) {
0285     resonanceWindow.resetCounter();
0286   }
0287   // Count events in background windows
0288   for (auto& backgroundWindow : backgroundWindow_) {
0289     backgroundWindow.resetCounter();
0290   }
0291
0292   // Now count the events in each window
0293   std::pair<lorentzVector, lorentzVector> muonPair;
0294   for (auto const& muonPair : muonPairs) {
0295     // Count events in resonance windows
0296     for (auto& resonanceWindow : resonanceWindow_) {
0297       resonanceWindow.count((muonPair.first + muonPair.second).mass(), weight);
0298     }
0299     // Count events in background windows
0300     for (auto& backgroundWindow : backgroundWindow_) {
0301       backgroundWindow.count((muonPair.first + muonPair.second).mass(), weight);
0302     }
0303   }
0304 }
0305
0306 void BackgroundHandler::consistencyCheck(const std::vector<int>& identifiers,
0307                                          const std::vector<double>& leftWindowBorders,
0308                                          const std::vector<double>& rightWindowBorders) const noexcept(false) {
0309   if (leftWindowBorders.size() != rightWindowBorders.size()) {
0310     throw cms::Exception("Configuration")
0311         << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = " << leftWindowBorders.size()
0312         << " != rightWindowBorders.size() = " << rightWindowBorders.size() << std::endl;
0313   }
0314   if (leftWindowBorders.size() != 3) {
0315     throw cms::Exception("Configuration")
0316         << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = rightWindowBorders.size() = "
0317         << leftWindowBorders.size() << " != 3" << std::endl;
0318   }
0319   if (identifiers.size() != 3) {
0320     throw cms::Exception("Configuration")
0321         << "BackgroundHandler::BackgroundHandler: identifiers must match the number of regions = 3" << std::endl;
0322   }
0323 }
0324
0325 #endif  // BackgroundHandler_cc