CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/MuonAnalysis/MomentumScaleCalibration/src/BackgroundHandler.cc

Go to the documentation of this file.
00001 #ifndef BackgroundHandler_cc
00002 #define BackgroundHandler_cc
00003 
00004 #include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundHandler.h"
00005 #include <algorithm>
00006 #include <TF1.h>
00007 #include <iostream>
00008 #include <boost/foreach.hpp>
00009 
00010 BackgroundHandler::BackgroundHandler( const std::vector<int> & identifiers,
00011                                       const std::vector<double> & leftWindowBorders,
00012                                       const std::vector<double> & rightWindowBorders,
00013                                       const double * ResMass,
00014                                       const double * massWindowHalfWidth )
00015 {
00016   // : leftWindowFactors_(leftWindowFactors),
00017   // rightWindowFactors_(rightWindowFactors)
00018 
00019   // Define the correspondence between regions and halfWidth to use
00020   // Defines also the function type to use (but they are checked to be consistent over a region)
00021   regToResHW_[0] = 0; // Region 0 use the one from Z
00022   regToResHW_[1] = 3; // Region 1 use the one from Upsilon1S
00023   regToResHW_[2] = 5; // Region 2 use the one from J/Psi
00024 
00025   // Define the correspondence between resonances and regions
00026   resToReg_[0] = 0; // Z
00027   resToReg_[1] = 1; // Upsilon3S
00028   resToReg_[2] = 1; // Upsilon2S
00029   resToReg_[3] = 1; // Upsilon1S
00030   resToReg_[4] = 2; // Psi2S
00031   resToReg_[5] = 2; // J/Psi
00032 
00033   // Throws cms::Exception("Configuration") in case the parameters are not what is expected
00034   consistencyCheck(identifiers, leftWindowBorders, rightWindowBorders);
00035 
00036   // Build the resonance windows
00037   for( unsigned int i=0; i<6; ++i ) {
00038     double mass = ResMass[i];
00039     double lowerLimit = mass - massWindowHalfWidth[i];
00040     double upperLimit = mass + massWindowHalfWidth[i];
00041     resonanceWindow_.push_back( MassWindow( mass, lowerLimit, upperLimit,
00042                                             std::vector<unsigned int>(1,i),
00043                                             backgroundFunctionService(identifiers[resToReg_[i]],
00044                                                                       lowerLimit,
00045                                                                       upperLimit) ) );
00046   }
00047 
00048   // Build the background windows
00049   // ----------------------------
00050   // Compute the mass center of each region
00051   double resMassForRegion[3];
00052   resMassForRegion[0] = ResMass[0];
00053   resMassForRegion[1] = (ResMass[1]+ResMass[2]+ResMass[3])/3;
00054   resMassForRegion[2] = (ResMass[4]+ResMass[5])/2;
00055 
00056   // Define which resonance is in which background window
00057   std::vector<std::vector<unsigned int> > indexes;
00058   indexes.push_back(std::vector<unsigned int>(1,0));
00059   indexes.push_back(std::vector<unsigned int>());
00060   for( int i=1; i<=3; ++i ) { indexes[1].push_back(i); }
00061   indexes.push_back(std::vector<unsigned int>());
00062   for( int i=4; i<=5; ++i ) { indexes[2].push_back(i); }
00063 
00064   unsigned int i=0;
00065   typedef std::vector<unsigned int> indexVec;
00066   BOOST_FOREACH(const indexVec & index, indexes) {
00067     //     double lowerLimit = resMassForRegion[i] - massWindowHalfWidth[regToResHW_[i]]*leftWindowFactors[i];
00068     //     double upperLimit = resMassForRegion[i] + massWindowHalfWidth[regToResHW_[i]]*rightWindowFactors[i];
00069     //     backgroundWindow_.push_back( MassWindow( resMassForRegion[i], lowerLimit, upperLimit, index,
00070     //                                              backgroundFunctionService(identifiers[i], lowerLimit, upperLimit ) ) );
00071     backgroundWindow_.push_back( MassWindow( resMassForRegion[i], leftWindowBorders[i], rightWindowBorders[i], index,
00072                                              backgroundFunctionService(identifiers[i], leftWindowBorders[i], rightWindowBorders[i] ) ) );
00073     ++i;
00074   }
00075   // Initialize the parNums to be used in the shifts of parval
00076   initializeParNums();
00077 }
00078 
00079 BackgroundHandler::~BackgroundHandler()
00080 {
00081 }
00082 
00083 void BackgroundHandler::initializeParNums()
00084 {
00085   // Initialize the parNums to be used in the shifts of parval
00086   parNumsRegions_[0] = 0;
00087   for( unsigned int i=1; i<backgroundWindow_.size() ; ++i ) {
00088     parNumsRegions_[i] = parNumsRegions_[i-1] + backgroundWindow_[i-1].backgroundFunction()->parNum();
00089   }
00090   parNumsResonances_[0] = parNumsRegions_[2]+backgroundWindow_[2].backgroundFunction()->parNum();
00091   for( unsigned int i=1; i<resonanceWindow_.size() ; ++i ) {
00092     parNumsResonances_[i] = parNumsResonances_[i-1] + resonanceWindow_[i-1].backgroundFunction()->parNum();
00093   }
00094 }
00095 
00096 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)
00097 {
00098   std::vector<double>::const_iterator parBgrIt = parBgr.begin();
00099   std::vector<int>::const_iterator parBgrOrderIt = parBgrOrder.begin();
00100   // Set the parameters for the regions only if this is not a rescaling
00101   for( int iReg = 0; iReg < 3; ++iReg ) {
00102     int shift = parNumsRegions_[iReg];
00103     backgroundWindow_[iReg].backgroundFunction()->setParameters( &(Start[shift]), &(Step[shift]), &(Mini[shift]),
00104         &(Maxi[shift]), &(ind[shift]), &(parname[shift]),
00105         parBgrIt+shift, parBgrOrderIt+shift, muonType );
00106   }
00107   for( int iRes = 0; iRes < 6; ++iRes ) {
00108     // parNumsResonances is already shifted for the regions parameters
00109     int shift = parNumsResonances_[iRes];
00110     resonanceWindow_[iRes].backgroundFunction()->setParameters( &(Start[shift]), &(Step[shift]), &(Mini[shift]),
00111         &(Maxi[shift]), &(ind[shift]), &(parname[shift]),
00112         parBgrIt+shift, parBgrOrderIt+shift, muonType );
00113   }
00114 }
00115 
00116 bool BackgroundHandler::unlockParameter(const std::vector<int> & resfind, const unsigned int ipar)
00117 {
00118   // parNumsRegions_ are shifted: [1] contains the number of parameters for 0 and so on.
00119   if( ipar < unsigned(parNumsRegions_[1]) && resfind[0] > 0 ) {
00120     return true;
00121   }
00122   if( ipar >= unsigned(parNumsRegions_[1]) && ipar < unsigned(parNumsRegions_[2]) && ( resfind[1] > 0 || resfind[2] > 0 || resfind[3] > 0 ) ) {
00123     return true;
00124   }
00125   // The first of parNumsResonances_ has the sum of parNums of the regions.
00126   if( ipar >= unsigned(parNumsRegions_[2]) && ipar < unsigned(parNumsResonances_[0]) && ( resfind[4] > 0 || resfind[5] > 0 ) ) {
00127     return true;
00128   }
00129   return false;
00130 }
00131 
00132 // std::pair<double, double> BackgroundHandler::windowFactors( const bool doBackgroundFit, const int ires )
00133 // {
00134 //   if( doBackgroundFit ) {
00135 //     // Fitting the background: use the regions
00136 //     return std::make_pair(leftWindowFactors_[resToReg_[ires]], rightWindowFactors_[resToReg_[ires]]);
00137 //   }
00138 //   else {
00139 //     // Not fitting the background: use the resonances
00140 //     return std::make_pair(1.,1.);
00141 //   }
00142 // }
00143 
00144 std::pair<double, double> BackgroundHandler::windowBorders( const bool doBackgroundFit, const int ires )
00145 {
00146   if( doBackgroundFit ) {
00147     // Fitting the background: use the regions
00148     return std::make_pair(backgroundWindow_[resToReg_[ires]].lowerBound(), backgroundWindow_[resToReg_[ires]].upperBound());
00149   }
00150   else {
00151     // Not fitting the background: use the resonances
00152     // return std::make_pair(1.,1.);
00153     return std::make_pair(resonanceWindow_[ires].lowerBound(), resonanceWindow_[ires].upperBound());
00154   }
00155 }
00156 
00157 double BackgroundHandler::resMass( const bool doBackgroundFit, const int ires )
00158 {
00159   if( doBackgroundFit ) {
00160     // Fitting the background: use the regions
00161     return backgroundWindow_[resToReg_[ires]].mass();
00162   }
00163   else {
00164     // Not fitting the background: use the resonances
00165     return resonanceWindow_[ires].mass();
00166   }
00167 }
00168 
00169 void BackgroundHandler::rescale( std::vector<double> & parBgr, const double * ResMass, const double * massWindowHalfWidth,
00170                                  const std::vector<std::pair<reco::Particle::LorentzVector,reco::Particle::LorentzVector> > & muonPairs,
00171                                  const double & weight )
00172 {
00173   countEventsInAllWindows(muonPairs, weight);
00174 
00175   // Loop on all regions and on all the resonances of each region and compute the background fraction
00176   // for each resonance window.
00177   unsigned int iRegion = 0;
00178   BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_)
00179   {
00180     // Iterator pointing to the first parameter of this background function in the full set of parameters
00181     std::vector<double>::const_iterator parBgrIt = (parBgr.begin()+parNumsRegions_[iRegion]);
00182     TF1 * backgroundFunctionForIntegral = backgroundWindow.backgroundFunction()->functionForIntegral(parBgrIt);
00183     // WARNING: this expects the background fraction parameter to be parBgr[0] for all the background functions.
00184     double kOld = *parBgrIt;
00185     double Nbw = backgroundWindow.events();
00186     double Ibw = backgroundFunctionForIntegral->Integral(backgroundWindow.lowerBound(),
00187                                                          backgroundWindow.upperBound());
00188 
00189     // index is the index of the resonance in the background window
00190     BOOST_FOREACH( unsigned int index, *(backgroundWindow.indexes()) )
00191     {
00192       // First set all parameters of the resonance window background function to those of the corresponding region
00193       for( int iPar = 0; iPar < resonanceWindow_[index].backgroundFunction()->parNum(); ++iPar ) {
00194         parBgr[parNumsResonances_[index]+iPar] = parBgr[parNumsRegions_[resToReg_[index]]+iPar];
00195       }
00196       // Estimated fraction of events in the resonance window
00197       double Irw = backgroundFunctionForIntegral->Integral(resonanceWindow_[index].lowerBound(),
00198                                                            resonanceWindow_[index].upperBound());
00199       double Nrw = resonanceWindow_[index].events();
00200 
00201       // Ibw is 1 (to avoid effects from approximation errors we set it to 1 and do not include it in the computation).
00202       if( Nrw != 0 ) parBgr[parNumsResonances_[index]] = kOld*Nbw/Nrw*Irw;
00203       else parBgr[parNumsResonances_[index]] = 0.;
00204 
00205       // Protect against fluctuations of number of events which could cause the fraction to go above 1.
00206       if( parBgr[parNumsResonances_[index]] > 1. ) parBgr[parNumsResonances_[index]] = 1.;
00207 
00208       double kNew = parBgr[parNumsResonances_[index]];
00209       std::cout << "For resonance = " << index << std::endl;
00210       std::cout << "backgroundWindow.lowerBound = " << backgroundWindow.lowerBound() << std::endl;
00211       std::cout << "backgroundWindow.upperBound = " << backgroundWindow.upperBound() << std::endl;
00212       std::cout << "parNumsResonances_["<<index<<"] = " << parNumsResonances_[index] << std::endl;
00213       std::cout << "Nbw = " << Nbw << ", Ibw = " << Ibw << std::endl;
00214       std::cout << "Nrw = " << Nrw << ", Irw = " << Irw << std::endl;
00215       std::cout << "k = " << kOld << ", k' = " << parBgr[parNumsResonances_[index]] << std::endl;
00216       std::cout << "background fraction in background window = Nbw*k = " << Nbw*kOld << std::endl;
00217       std::cout << "background fraction in resonance window = Nrw*k' = " << Nrw*kNew << std::endl;
00218     }
00219     ++iRegion;
00220     delete backgroundFunctionForIntegral;
00221   }
00222 }
00223 
00224 std::pair<double, double> BackgroundHandler::backgroundFunction( const bool doBackgroundFit,
00225                                                                  const double * parval, const int resTotNum, const int ires,
00226                                                                  const bool * resConsidered, const double * ResMass, const double ResHalfWidth[],
00227                                                                  const int MuonType, const double & mass, const double & resEta )
00228 {
00229   if( doBackgroundFit ) {
00230     // Return the values for the region
00231     int iReg = resToReg_[ires];
00232     return std::make_pair( parval[parNumsRegions_[iReg]] * backgroundWindow_[iReg].backgroundFunction()->fracVsEta(&(parval[parNumsRegions_[iReg]]), resEta),
00233                            (*(backgroundWindow_[iReg].backgroundFunction()))( &(parval[parNumsRegions_[iReg]]), mass, resEta ) );
00234   }
00235   // Return the values for the resonance
00236   return std::make_pair( parval[parNumsResonances_[ires]] * resonanceWindow_[ires].backgroundFunction()->fracVsEta(&(parval[parNumsResonances_[ires]]), resEta),
00237                          (*(resonanceWindow_[ires].backgroundFunction()))( &(parval[parNumsResonances_[ires]]), mass, resEta ) );
00238 }
00239 
00240 void BackgroundHandler::countEventsInAllWindows(const std::vector<std::pair<reco::Particle::LorentzVector,reco::Particle::LorentzVector> > & muonPairs,
00241                                                 const double & weight)
00242 {
00243   // First reset all the counters
00244   BOOST_FOREACH(MassWindow & resonanceWindow, resonanceWindow_) {
00245     resonanceWindow.resetCounter();
00246   }
00247   // Count events in background windows
00248   BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_) {
00249     backgroundWindow.resetCounter();
00250   }
00251 
00252   // Now count the events in each window
00253   std::pair<lorentzVector,lorentzVector> muonPair;
00254   BOOST_FOREACH(muonPair, muonPairs) {
00255     // Count events in resonance windows
00256     BOOST_FOREACH(MassWindow & resonanceWindow, resonanceWindow_) {
00257       resonanceWindow.count((muonPair.first + muonPair.second).mass(), weight);
00258     }
00259     // Count events in background windows
00260     BOOST_FOREACH(MassWindow & backgroundWindow, backgroundWindow_) {
00261       backgroundWindow.count((muonPair.first + muonPair.second).mass(), weight);
00262     }
00263   }
00264 }
00265 
00266 void BackgroundHandler::consistencyCheck(const std::vector<int> & identifiers,
00267                                          const std::vector<double> & leftWindowBorders,
00268                                          const std::vector<double> & rightWindowBorders) const throw(cms::Exception)
00269 {
00270   if( leftWindowBorders.size() != rightWindowBorders.size() ) {
00271     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = " << leftWindowBorders.size()
00272                                           << " != rightWindowBorders.size() = " << rightWindowBorders.size() << std::endl;
00273   }
00274   if( leftWindowBorders.size() != 3 ) {
00275     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: leftWindowBorders.size() = rightWindowBorders.size() = "
00276                                           << leftWindowBorders.size() << " != 3" << std::endl;
00277   }
00278   if( identifiers.size() != 3 ) {
00279     throw cms::Exception("Configuration") << "BackgroundHandler::BackgroundHandler: identifiers must match the number of regions = 3" << std::endl;
00280   }
00281 }
00282 
00283 #endif // BackgroundHandler_cc