Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:29:33

0001 /** \class AlcaBeamSpotManager
0002  *  No description available.
0003  *
0004  *  \author L. Uplegger F. Yumiceva - Fermilab
0005  */
0006 
0007 #include "Calibration/TkAlCaRecoProducers/interface/AlcaBeamSpotManager.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include <climits>
0011 #include <cmath>
0012 #include <vector>
0013 
0014 using namespace edm;
0015 using namespace reco;
0016 using namespace std;
0017 
0018 //--------------------------------------------------------------------------------------------------
0019 AlcaBeamSpotManager::AlcaBeamSpotManager(void) {}
0020 
0021 //--------------------------------------------------------------------------------------------------
0022 AlcaBeamSpotManager::AlcaBeamSpotManager(const ParameterSet &iConfig, edm::ConsumesCollector &&iC)
0023     : beamSpotOutputBase_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters")
0024                               .getUntrackedParameter<std::string>("BeamSpotOutputBase")),
0025       beamSpotModuleName_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters")
0026                               .getUntrackedParameter<std::string>("BeamSpotModuleName")),
0027       beamSpotLabel_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters")
0028                          .getUntrackedParameter<std::string>("BeamSpotLabel")),
0029       sigmaZCut_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters")
0030                      .getUntrackedParameter<double>("SigmaZCut")) {
0031   edm::InputTag beamSpotTag_(beamSpotModuleName_, beamSpotLabel_);
0032   beamSpotToken_ = iC.consumes<reco::BeamSpot, edm::InLumi>(beamSpotTag_);
0033   LogInfo("AlcaBeamSpotManager") << "Output base: " << beamSpotOutputBase_ << std::endl;
0034   reset();
0035 }
0036 
0037 //--------------------------------------------------------------------------------------------------
0038 AlcaBeamSpotManager::~AlcaBeamSpotManager(void) {}
0039 
0040 //--------------------------------------------------------------------------------------------------
0041 void AlcaBeamSpotManager::reset(void) { beamSpotMap_.clear(); }
0042 //--------------------------------------------------------------------------------------------------
0043 void AlcaBeamSpotManager::readLumi(const LuminosityBlock &iLumi) {
0044   Handle<BeamSpot> beamSpotHandle;
0045   iLumi.getByToken(beamSpotToken_, beamSpotHandle);
0046 
0047   if (beamSpotHandle.isValid()) {  // check the product
0048     std::pair<edm::Timestamp, reco::BeamSpot> time_bs(iLumi.beginTime(), *beamSpotHandle);
0049     beamSpotMap_[iLumi.luminosityBlock()] = time_bs;
0050     const BeamSpot *aBeamSpot = &beamSpotMap_[iLumi.luminosityBlock()].second;
0051     aBeamSpot = beamSpotHandle.product();
0052     LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
0053     LogInfo("AlcaBeamSpotManager") << *aBeamSpot << std::endl;
0054   } else {
0055     LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
0056     LogInfo("AlcaBeamSpotManager") << "   BS is not valid!" << std::endl;
0057   }
0058 }
0059 
0060 //--------------------------------------------------------------------------------------------------
0061 void AlcaBeamSpotManager::createWeightedPayloads(void) {
0062   vector<bsMap_iterator> listToErase;
0063   for (bsMap_iterator it = beamSpotMap_.begin(); it != beamSpotMap_.end(); it++) {
0064     if (it->second.second.type() != BeamSpot::Tracker || it->second.second.sigmaZ() < sigmaZCut_) {
0065       listToErase.push_back(it);
0066     }
0067   }
0068   for (vector<bsMap_iterator>::iterator it = listToErase.begin(); it != listToErase.end(); it++) {
0069     beamSpotMap_.erase(*it);
0070   }
0071   if (beamSpotMap_.size() <= 1) {
0072     return;
0073   }
0074   // Return only if lumibased since the collapsing alghorithm requires the next
0075   // and next to next lumi sections
0076   else if (beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased") {
0077     return;
0078   }
0079   if (beamSpotOutputBase_ == "lumibased") {
0080     //    bsMap_iterator referenceBS = beamSpotMap_.begin();
0081     bsMap_iterator firstBS = beamSpotMap_.begin();
0082     //    bsMap_iterator lastBS      = beamSpotMap_.begin();
0083     bsMap_iterator currentBS = beamSpotMap_.begin();
0084     bsMap_iterator nextBS = ++beamSpotMap_.begin();
0085     bsMap_iterator nextNextBS = ++(++(beamSpotMap_.begin()));
0086 
0087     reco::BeamSpot currentBSObj;
0088     reco::BeamSpot nextBSObj;
0089 
0090     map<LuminosityBlockNumber_t, std::pair<edm::Timestamp, BeamSpot>> tmpBeamSpotMap_;
0091     bool docreate = true;
0092     bool endOfRun = false;  // Added
0093     bool docheck = true;    // Added
0094     bool foundShift = false;
0095     long countlumi = 0;   // Added
0096     string tmprun = "";   // Added
0097     long maxNlumis = 20;  // Added
0098     //    if weighted:
0099     //        maxNlumis = 999999999
0100 
0101     unsigned int iteration = 0;
0102     // while(nextNextBS!=beamSpotMap_.end()){
0103     while (nextBS != beamSpotMap_.end()) {
0104       LogInfo("AlcaBeamSpotManager") << "Iteration: " << iteration << " size: " << beamSpotMap_.size() << "\n"
0105                                      << "Lumi: " << currentBS->first << "\n"
0106                                      << currentBS->second.second << "\n"
0107                                      << nextBS->first << "\n"
0108                                      << nextBS->second.second << endl;
0109       if (nextNextBS != beamSpotMap_.end())
0110         LogInfo("AlcaBeamSpotManager") << nextNextBS->first << "\n" << nextNextBS->second.second << endl;
0111 
0112       if (docreate) {
0113         firstBS = currentBS;
0114         docreate = false;  // Added
0115       }
0116       // if(iteration >= beamSpotMap_.size()-3){
0117       if (iteration >= beamSpotMap_.size() - 2) {
0118         LogInfo("AlcaBeamSpotManager") << "Reached lumi " << currentBS->first
0119                                        << " now close payload because end of data has been reached.";
0120         docreate = true;
0121         endOfRun = true;
0122       }
0123       // check we run over the same run
0124       //      if (ibeam->first.Run() != inextbeam->first.Run()){
0125       //        LogInfo("AlcaBeamSpotManager")
0126       //          << "close payload because end of run.";
0127       //        docreate = true;
0128       //      }
0129       // check maximum lumi counts
0130       if (countlumi == maxNlumis - 1) {
0131         LogInfo("AlcaBeamSpotManager") << "close payload because maximum lumi "
0132                                           "sections accumulated within run ";
0133         docreate = true;
0134         countlumi = 0;
0135       }
0136       //      if weighted:
0137       //          docheck = False
0138       // check offsets
0139       if (docheck) {
0140         foundShift = false;
0141         LogInfo("AlcaBeamSpotManager") << "Checking checking!" << endl;
0142         float limit = 0;
0143         pair<float, float> adelta1;
0144         pair<float, float> adelta2;
0145         pair<float, float> adelta;
0146         pair<float, float> adelta1dxdz;
0147         pair<float, float> adelta2dxdz;
0148         pair<float, float> adelta1dydz;
0149         pair<float, float> adelta2dydz;
0150         pair<float, float> adelta1widthX;
0151         pair<float, float> adelta2widthX;
0152         pair<float, float> adelta1widthY;
0153         pair<float, float> adelta2widthY;
0154         pair<float, float> adelta1z0;
0155         pair<float, float> adelta1sigmaZ;
0156 
0157         // define minimum limit
0158         float min_limit = 0.0025;
0159 
0160         // limit for x and y
0161         limit = currentBS->second.second.BeamWidthX() / 2.;
0162         if (limit < min_limit) {
0163           limit = min_limit;
0164         }
0165 
0166         currentBSObj = currentBS->second.second;
0167         nextBSObj = nextBS->second.second;
0168 
0169         // check movements in X
0170         adelta1 = delta(currentBSObj.x0(), currentBSObj.x0Error(), nextBSObj.x0(), nextBSObj.x0Error());
0171         adelta2 = pair<float, float>(0., 1.e9);
0172         if (nextNextBS->second.second.type() != -1) {
0173           adelta2 = delta(nextBS->second.second.x0(),
0174                           nextBSObj.x0Error(),
0175                           nextNextBS->second.second.x0(),
0176                           nextNextBS->second.second.x0Error());
0177         }
0178         bool deltaX = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
0179         if (iteration < beamSpotMap_.size() - 2) {
0180           if (!deltaX && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
0181             LogInfo("AlcaBeamSpotManager")
0182                 << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
0183             deltaX = true;
0184           } else if (deltaX && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
0185                      fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
0186             LogInfo("AlcaBeamSpotManager") << " negative, " << adelta1.first / adelta2.first << endl;
0187             deltaX = false;
0188           }
0189         }
0190 
0191         // calculating all deltas
0192         adelta1dxdz = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
0193         adelta2dxdz = pair<float, float>(0., 1.e9);
0194         adelta1dydz = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
0195         adelta2dydz = pair<float, float>(0., 1.e9);
0196         adelta1widthX = delta(currentBSObj.BeamWidthX(),
0197                               currentBSObj.BeamWidthXError(),
0198                               nextBSObj.BeamWidthX(),
0199                               nextBSObj.BeamWidthXError());
0200         adelta2widthX = pair<float, float>(0., 1.e9);
0201         adelta1widthY = delta(currentBSObj.BeamWidthY(),
0202                               currentBSObj.BeamWidthYError(),
0203                               nextBSObj.BeamWidthY(),
0204                               nextBSObj.BeamWidthYError());
0205         adelta2widthY = pair<float, float>(0., 1.e9);
0206         adelta1z0 = delta(currentBSObj.z0(), currentBSObj.z0Error(), nextBSObj.z0(), nextBSObj.z0Error());
0207         adelta1sigmaZ =
0208             delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
0209 
0210         // check movements in Y
0211         adelta1 = delta(currentBSObj.y0(), currentBSObj.y0Error(), nextBSObj.y0(), nextBSObj.y0Error());
0212         adelta2 = pair<float, float>(0., 1.e9);
0213         if (nextNextBS->second.second.type() != BeamSpot::Unknown) {
0214           adelta2 = delta(
0215               nextBSObj.y0(), nextBSObj.y0Error(), nextNextBS->second.second.y0(), nextNextBS->second.second.y0Error());
0216           adelta2dxdz = delta(nextBSObj.dxdz(),
0217                               nextBSObj.dxdzError(),
0218                               nextNextBS->second.second.dxdz(),
0219                               nextNextBS->second.second.dxdzError());
0220           adelta2dydz = delta(nextBSObj.dydz(),
0221                               nextBSObj.dydzError(),
0222                               nextNextBS->second.second.dydz(),
0223                               nextNextBS->second.second.dydzError());
0224           adelta2widthX = delta(nextBSObj.BeamWidthX(),
0225                                 nextBSObj.BeamWidthXError(),
0226                                 nextNextBS->second.second.BeamWidthX(),
0227                                 nextNextBS->second.second.BeamWidthXError());
0228           adelta2widthY = delta(nextBSObj.BeamWidthY(),
0229                                 nextBSObj.BeamWidthYError(),
0230                                 nextNextBS->second.second.BeamWidthY(),
0231                                 nextNextBS->second.second.BeamWidthYError());
0232         }
0233         bool deltaY = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
0234         if (iteration < beamSpotMap_.size() - 2) {
0235           if (!deltaY && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
0236             LogInfo("AlcaBeamSpotManager")
0237                 << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
0238             deltaY = true;
0239           } else if (deltaY && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
0240                      fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
0241             LogInfo("AlcaBeamSpotManager") << " negative, " << adelta1.first / adelta2.first << endl;
0242             deltaY = false;
0243           }
0244         }
0245 
0246         limit = currentBSObj.sigmaZ() / 2.;
0247         bool deltaZ =
0248             (deltaSig(adelta1z0.first, adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >= limit) ? true : false;
0249         adelta =
0250             delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
0251         bool deltasigmaZ = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0252         bool deltadxdz = false;
0253         bool deltadydz = false;
0254         bool deltawidthX = false;
0255         bool deltawidthY = false;
0256 
0257         if (iteration < beamSpotMap_.size() - 2) {
0258           adelta = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
0259           deltadxdz = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0260           if (deltadxdz && (adelta1dxdz.first * adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 &&
0261               fabs(adelta1dxdz.first / adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first / adelta2dxdz.first) < 3) {
0262             deltadxdz = false;
0263           }
0264 
0265           adelta = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
0266           deltadydz = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0267           if (deltadydz && (adelta1dydz.first * adelta2dydz.first) < 0 && adelta2dydz.first != 0 &&
0268               fabs(adelta1dydz.first / adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first / adelta2dydz.first) < 3) {
0269             deltadydz = false;
0270           }
0271 
0272           adelta = delta(currentBSObj.BeamWidthX(),
0273                          currentBSObj.BeamWidthXError(),
0274                          nextBSObj.BeamWidthX(),
0275                          nextBSObj.BeamWidthXError());
0276           deltawidthX = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0277           if (deltawidthX && (adelta1widthX.first * adelta2widthX.first) < 0 && adelta2widthX.first != 0 &&
0278               fabs(adelta1widthX.first / adelta2widthX.first) > 0.33 &&
0279               fabs(adelta1widthX.first / adelta2widthX.first) < 3) {
0280             deltawidthX = false;
0281           }
0282 
0283           adelta = delta(currentBSObj.BeamWidthY(),
0284                          currentBSObj.BeamWidthYError(),
0285                          nextBSObj.BeamWidthY(),
0286                          nextBSObj.BeamWidthYError());
0287           deltawidthY = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0288           if (deltawidthY && (adelta1widthY.first * adelta2widthY.first) < 0 && adelta2widthY.first != 0 &&
0289               fabs(adelta1widthY.first / adelta2widthY.first) > 0.33 &&
0290               fabs(adelta1widthY.first / adelta2widthY.first) < 3) {
0291             deltawidthY = false;
0292           }
0293         }
0294         if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY) {
0295           docreate = true;
0296           foundShift = true;
0297           LogInfo("AlcaBeamSpotManager") << "close payload because of movement in"
0298                                          << " X=" << deltaX << ", Y=" << deltaY << ", Z=" << deltaZ
0299                                          << ", sigmaZ=" << deltasigmaZ << ", dxdz=" << deltadxdz
0300                                          << ", dydz=" << deltadydz << ", widthX=" << deltawidthX
0301                                          << ", widthY=" << deltawidthY << endl;
0302         }
0303       }
0304       if (docreate) {
0305         std::pair<edm::Timestamp, reco::BeamSpot> thepair;
0306         if (foundShift) {
0307           thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
0308           tmpBeamSpotMap_[firstBS->first] = thepair;
0309           if (endOfRun) {
0310             // if we're here, then we need to found a shift in the last LS
0311             // We already created a new IOV, now create one just for the last LS
0312             thepair = std::make_pair(nextBS->second.first, nextBSObj);
0313             tmpBeamSpotMap_[nextBS->first] = thepair;
0314           }
0315         } else if (!foundShift && !endOfRun) {  // maxLS reached
0316           thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
0317           tmpBeamSpotMap_[firstBS->first] = thepair;
0318         } else {  // end of run with no shift detectred in last LS
0319           thepair = std::make_pair(currentBS->second.first, weight(firstBS, beamSpotMap_.end()));
0320           tmpBeamSpotMap_[firstBS->first] = thepair;
0321         }
0322         firstBS = nextBS;
0323         countlumi = 0;
0324       }
0325       // tmprun = currentBS->second.Run
0326       // increase the counter by one only if the IOV hasn't been closed
0327       if (!docreate)
0328         ++countlumi;
0329 
0330       currentBS = nextBS;
0331       nextBS = nextNextBS;
0332       nextNextBS++;
0333       ++iteration;
0334     }
0335     beamSpotMap_.clear();
0336     beamSpotMap_ = tmpBeamSpotMap_;
0337   } else if (beamSpotOutputBase_ == "runbased") {
0338     LuminosityBlockNumber_t firstLumi = beamSpotMap_.begin()->first;
0339     std::pair<edm::Timestamp, reco::BeamSpot> thepair(beamSpotMap_.begin()->second.first,
0340                                                       weight(beamSpotMap_.begin(), beamSpotMap_.end()));
0341     beamSpotMap_.clear();
0342     beamSpotMap_[firstLumi] = thepair;
0343   } else {
0344     LogInfo("AlcaBeamSpotManager") << "Unrecognized BeamSpotOutputBase parameter: " << beamSpotOutputBase_ << endl;
0345   }
0346 }
0347 
0348 //--------------------------------------------------------------------------------------------------
0349 BeamSpot AlcaBeamSpotManager::weight(const bsMap_iterator &begin, const bsMap_iterator &end) {
0350   double x, xError = 0;
0351   double y, yError = 0;
0352   double z, zError = 0;
0353   double sigmaZ, sigmaZError = 0;
0354   double dxdz, dxdzError = 0;
0355   double dydz, dydzError = 0;
0356   double widthX, widthXError = 0;
0357   double widthY, widthYError = 0;
0358   LogInfo("AlcaBeamSpotManager") << "Weighted BeamSpot will span lumi " << begin->first << " to " << end->first << endl;
0359 
0360   BeamSpot::BeamType type = BeamSpot::Unknown;
0361   for (bsMap_iterator it = begin; it != end; it++) {
0362     weight(x, xError, it->second.second.x0(), it->second.second.x0Error());
0363     weight(y, yError, it->second.second.y0(), it->second.second.y0Error());
0364     weight(z, zError, it->second.second.z0(), it->second.second.z0Error());
0365     weight(sigmaZ, sigmaZError, it->second.second.sigmaZ(), it->second.second.sigmaZ0Error());
0366     weight(dxdz, dxdzError, it->second.second.dxdz(), it->second.second.dxdzError());
0367     weight(dydz, dydzError, it->second.second.dydz(), it->second.second.dydzError());
0368     weight(widthX, widthXError, it->second.second.BeamWidthX(), it->second.second.BeamWidthXError());
0369     weight(widthY, widthYError, it->second.second.BeamWidthY(), it->second.second.BeamWidthYError());
0370     if (it->second.second.type() == BeamSpot::Tracker) {
0371       type = BeamSpot::Tracker;
0372     }
0373   }
0374   BeamSpot::Point bsPosition(x, y, z);
0375   BeamSpot::CovarianceMatrix error;
0376   error(0, 0) = xError * xError;
0377   error(1, 1) = yError * yError;
0378   error(2, 2) = zError * zError;
0379   error(3, 3) = sigmaZError * sigmaZError;
0380   error(4, 4) = dxdzError * dxdzError;
0381   error(5, 5) = dydzError * dydzError;
0382   error(6, 6) = widthXError * widthXError;
0383   BeamSpot weightedBeamSpot(bsPosition, sigmaZ, dxdz, dydz, widthX, error, type);
0384   weightedBeamSpot.setBeamWidthY(widthY);
0385   LogInfo("AlcaBeamSpotManager") << "Weighted BeamSpot will be:" << '\n' << weightedBeamSpot << endl;
0386   return weightedBeamSpot;
0387 }
0388 
0389 //--------------------------------------------------------------------------------------------------
0390 void AlcaBeamSpotManager::weight(double &mean, double &meanError, const double &val, const double &valError) {
0391   double tmpError = 0;
0392   if (meanError < 1e-8) {
0393     tmpError = 1 / (valError * valError);
0394     mean = val * tmpError;
0395   } else {
0396     tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
0397     mean = mean / (meanError * meanError) + val / (valError * valError);
0398   }
0399   mean = mean / tmpError;
0400   meanError = sqrt(1 / tmpError);
0401 }
0402 
0403 //--------------------------------------------------------------------------------------------------
0404 pair<float, float> AlcaBeamSpotManager::delta(const float &x,
0405                                               const float &xError,
0406                                               const float &nextX,
0407                                               const float &nextXError) {
0408   return pair<float, float>(x - nextX, sqrt(pow(xError, 2) + pow(nextXError, 2)));
0409 }
0410 
0411 //--------------------------------------------------------------------------------------------------
0412 float AlcaBeamSpotManager::deltaSig(const float &num, const float &den) {
0413   if (den != 0) {
0414     return fabs(num / den);
0415   } else {
0416     return float(LONG_MAX);
0417   }
0418 }