File indexing completed on 2025-03-27 03:36:13
0001
0002
0003
0004
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()) {
0048 std::pair<edm::Timestamp, reco::BeamSpot> time_bs(iLumi.beginTime(), *beamSpotHandle);
0049 beamSpotMap_[iLumi.luminosityBlock()] = time_bs;
0050 const BeamSpot *aBeamSpot = beamSpotHandle.product();
0051 LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
0052 LogInfo("AlcaBeamSpotManager") << *aBeamSpot << std::endl;
0053 } else {
0054 LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
0055 LogInfo("AlcaBeamSpotManager") << " BS is not valid!" << std::endl;
0056 }
0057 }
0058
0059
0060 void AlcaBeamSpotManager::createWeightedPayloads(void) {
0061 vector<bsMap_iterator> listToErase;
0062 for (bsMap_iterator it = beamSpotMap_.begin(); it != beamSpotMap_.end(); it++) {
0063 if (it->second.second.type() != BeamSpot::Tracker || it->second.second.sigmaZ() < sigmaZCut_) {
0064 listToErase.push_back(it);
0065 }
0066 }
0067 for (vector<bsMap_iterator>::iterator it = listToErase.begin(); it != listToErase.end(); it++) {
0068 beamSpotMap_.erase(*it);
0069 }
0070 if (beamSpotMap_.size() <= 1) {
0071 return;
0072 }
0073
0074
0075 else if (beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased") {
0076 return;
0077 }
0078 if (beamSpotOutputBase_ == "lumibased") {
0079
0080 bsMap_iterator firstBS = beamSpotMap_.begin();
0081
0082 bsMap_iterator currentBS = beamSpotMap_.begin();
0083 bsMap_iterator nextBS = ++beamSpotMap_.begin();
0084 bsMap_iterator nextNextBS = ++(++(beamSpotMap_.begin()));
0085
0086 reco::BeamSpot currentBSObj;
0087 reco::BeamSpot nextBSObj;
0088
0089 map<LuminosityBlockNumber_t, std::pair<edm::Timestamp, BeamSpot>> tmpBeamSpotMap_;
0090 bool docreate = true;
0091 bool endOfRun = false;
0092 bool docheck = true;
0093 bool foundShift = false;
0094 long countlumi = 0;
0095 string tmprun = "";
0096 long maxNlumis = 20;
0097
0098
0099
0100 unsigned int iteration = 0;
0101
0102 while (nextBS != beamSpotMap_.end()) {
0103 LogInfo("AlcaBeamSpotManager") << "Iteration: " << iteration << " size: " << beamSpotMap_.size() << "\n"
0104 << "Lumi: " << currentBS->first << "\n"
0105 << currentBS->second.second << "\n"
0106 << nextBS->first << "\n"
0107 << nextBS->second.second << endl;
0108 if (nextNextBS != beamSpotMap_.end())
0109 LogInfo("AlcaBeamSpotManager") << nextNextBS->first << "\n" << nextNextBS->second.second << endl;
0110
0111 if (docreate) {
0112 firstBS = currentBS;
0113 docreate = false;
0114 }
0115
0116 if (iteration >= beamSpotMap_.size() - 2) {
0117 LogInfo("AlcaBeamSpotManager") << "Reached lumi " << currentBS->first
0118 << " now close payload because end of data has been reached.";
0119 docreate = true;
0120 endOfRun = true;
0121 }
0122
0123
0124
0125
0126
0127
0128
0129 if (countlumi == maxNlumis - 1) {
0130 LogInfo("AlcaBeamSpotManager") << "close payload because maximum lumi "
0131 "sections accumulated within run ";
0132 docreate = true;
0133 countlumi = 0;
0134 }
0135
0136
0137
0138 if (docheck) {
0139 foundShift = false;
0140 LogInfo("AlcaBeamSpotManager") << "Checking checking!" << endl;
0141 float limit = 0;
0142 pair<float, float> adelta1;
0143 pair<float, float> adelta2;
0144 pair<float, float> adelta;
0145 pair<float, float> adelta1dxdz;
0146 pair<float, float> adelta2dxdz;
0147 pair<float, float> adelta1dydz;
0148 pair<float, float> adelta2dydz;
0149 pair<float, float> adelta1widthX;
0150 pair<float, float> adelta2widthX;
0151 pair<float, float> adelta1widthY;
0152 pair<float, float> adelta2widthY;
0153 pair<float, float> adelta1z0;
0154 pair<float, float> adelta1sigmaZ;
0155
0156
0157 float min_limit = 0.0025;
0158
0159
0160 limit = currentBS->second.second.BeamWidthX() / 2.;
0161 if (limit < min_limit) {
0162 limit = min_limit;
0163 }
0164
0165 currentBSObj = currentBS->second.second;
0166 nextBSObj = nextBS->second.second;
0167
0168
0169 adelta1 = delta(currentBSObj.x0(), currentBSObj.x0Error(), nextBSObj.x0(), nextBSObj.x0Error());
0170 adelta2 = pair<float, float>(0., 1.e9);
0171 if (nextNextBS->second.second.type() != -1) {
0172 adelta2 = delta(nextBS->second.second.x0(),
0173 nextBSObj.x0Error(),
0174 nextNextBS->second.second.x0(),
0175 nextNextBS->second.second.x0Error());
0176 }
0177 bool deltaX = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
0178 if (iteration < beamSpotMap_.size() - 2) {
0179 if (!deltaX && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
0180 LogInfo("AlcaBeamSpotManager")
0181 << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
0182 deltaX = true;
0183 } else if (deltaX && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
0184 fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
0185 LogInfo("AlcaBeamSpotManager") << " negative, " << adelta1.first / adelta2.first << endl;
0186 deltaX = false;
0187 }
0188 }
0189
0190
0191 adelta1dxdz = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
0192 adelta2dxdz = pair<float, float>(0., 1.e9);
0193 adelta1dydz = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
0194 adelta2dydz = pair<float, float>(0., 1.e9);
0195 adelta1widthX = delta(currentBSObj.BeamWidthX(),
0196 currentBSObj.BeamWidthXError(),
0197 nextBSObj.BeamWidthX(),
0198 nextBSObj.BeamWidthXError());
0199 adelta2widthX = pair<float, float>(0., 1.e9);
0200 adelta1widthY = delta(currentBSObj.BeamWidthY(),
0201 currentBSObj.BeamWidthYError(),
0202 nextBSObj.BeamWidthY(),
0203 nextBSObj.BeamWidthYError());
0204 adelta2widthY = pair<float, float>(0., 1.e9);
0205 adelta1z0 = delta(currentBSObj.z0(), currentBSObj.z0Error(), nextBSObj.z0(), nextBSObj.z0Error());
0206 adelta1sigmaZ =
0207 delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
0208
0209
0210 adelta1 = delta(currentBSObj.y0(), currentBSObj.y0Error(), nextBSObj.y0(), nextBSObj.y0Error());
0211 adelta2 = pair<float, float>(0., 1.e9);
0212 if (nextNextBS->second.second.type() != BeamSpot::Unknown) {
0213 adelta2 = delta(
0214 nextBSObj.y0(), nextBSObj.y0Error(), nextNextBS->second.second.y0(), nextNextBS->second.second.y0Error());
0215 adelta2dxdz = delta(nextBSObj.dxdz(),
0216 nextBSObj.dxdzError(),
0217 nextNextBS->second.second.dxdz(),
0218 nextNextBS->second.second.dxdzError());
0219 adelta2dydz = delta(nextBSObj.dydz(),
0220 nextBSObj.dydzError(),
0221 nextNextBS->second.second.dydz(),
0222 nextNextBS->second.second.dydzError());
0223 adelta2widthX = delta(nextBSObj.BeamWidthX(),
0224 nextBSObj.BeamWidthXError(),
0225 nextNextBS->second.second.BeamWidthX(),
0226 nextNextBS->second.second.BeamWidthXError());
0227 adelta2widthY = delta(nextBSObj.BeamWidthY(),
0228 nextBSObj.BeamWidthYError(),
0229 nextNextBS->second.second.BeamWidthY(),
0230 nextNextBS->second.second.BeamWidthYError());
0231 }
0232 bool deltaY = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
0233 if (iteration < beamSpotMap_.size() - 2) {
0234 if (!deltaY && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
0235 LogInfo("AlcaBeamSpotManager")
0236 << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
0237 deltaY = true;
0238 } else if (deltaY && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
0239 fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
0240 LogInfo("AlcaBeamSpotManager") << " negative, " << adelta1.first / adelta2.first << endl;
0241 deltaY = false;
0242 }
0243 }
0244
0245 limit = currentBSObj.sigmaZ() / 2.;
0246 bool deltaZ =
0247 (deltaSig(adelta1z0.first, adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >= limit) ? true : false;
0248 adelta =
0249 delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
0250 bool deltasigmaZ = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0251 bool deltadxdz = false;
0252 bool deltadydz = false;
0253 bool deltawidthX = false;
0254 bool deltawidthY = false;
0255
0256 if (iteration < beamSpotMap_.size() - 2) {
0257 adelta = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
0258 deltadxdz = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0259 if (deltadxdz && (adelta1dxdz.first * adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 &&
0260 fabs(adelta1dxdz.first / adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first / adelta2dxdz.first) < 3) {
0261 deltadxdz = false;
0262 }
0263
0264 adelta = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
0265 deltadydz = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0266 if (deltadydz && (adelta1dydz.first * adelta2dydz.first) < 0 && adelta2dydz.first != 0 &&
0267 fabs(adelta1dydz.first / adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first / adelta2dydz.first) < 3) {
0268 deltadydz = false;
0269 }
0270
0271 adelta = delta(currentBSObj.BeamWidthX(),
0272 currentBSObj.BeamWidthXError(),
0273 nextBSObj.BeamWidthX(),
0274 nextBSObj.BeamWidthXError());
0275 deltawidthX = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0276 if (deltawidthX && (adelta1widthX.first * adelta2widthX.first) < 0 && adelta2widthX.first != 0 &&
0277 fabs(adelta1widthX.first / adelta2widthX.first) > 0.33 &&
0278 fabs(adelta1widthX.first / adelta2widthX.first) < 3) {
0279 deltawidthX = false;
0280 }
0281
0282 adelta = delta(currentBSObj.BeamWidthY(),
0283 currentBSObj.BeamWidthYError(),
0284 nextBSObj.BeamWidthY(),
0285 nextBSObj.BeamWidthYError());
0286 deltawidthY = (deltaSig(adelta.first, adelta.second) > 5.0) ? true : false;
0287 if (deltawidthY && (adelta1widthY.first * adelta2widthY.first) < 0 && adelta2widthY.first != 0 &&
0288 fabs(adelta1widthY.first / adelta2widthY.first) > 0.33 &&
0289 fabs(adelta1widthY.first / adelta2widthY.first) < 3) {
0290 deltawidthY = false;
0291 }
0292 }
0293 if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY) {
0294 docreate = true;
0295 foundShift = true;
0296 LogInfo("AlcaBeamSpotManager") << "close payload because of movement in"
0297 << " X=" << deltaX << ", Y=" << deltaY << ", Z=" << deltaZ
0298 << ", sigmaZ=" << deltasigmaZ << ", dxdz=" << deltadxdz
0299 << ", dydz=" << deltadydz << ", widthX=" << deltawidthX
0300 << ", widthY=" << deltawidthY << endl;
0301 }
0302 }
0303 if (docreate) {
0304 std::pair<edm::Timestamp, reco::BeamSpot> thepair;
0305 if (foundShift) {
0306 thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
0307 tmpBeamSpotMap_[firstBS->first] = thepair;
0308 if (endOfRun) {
0309
0310
0311 thepair = std::make_pair(nextBS->second.first, nextBSObj);
0312 tmpBeamSpotMap_[nextBS->first] = thepair;
0313 }
0314 } else if (!foundShift && !endOfRun) {
0315 thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
0316 tmpBeamSpotMap_[firstBS->first] = thepair;
0317 } else {
0318 thepair = std::make_pair(currentBS->second.first, weight(firstBS, beamSpotMap_.end()));
0319 tmpBeamSpotMap_[firstBS->first] = thepair;
0320 }
0321 firstBS = nextBS;
0322 countlumi = 0;
0323 }
0324
0325
0326 if (!docreate)
0327 ++countlumi;
0328
0329 currentBS = nextBS;
0330 nextBS = nextNextBS;
0331 nextNextBS++;
0332 ++iteration;
0333 }
0334 beamSpotMap_.clear();
0335 beamSpotMap_ = tmpBeamSpotMap_;
0336 } else if (beamSpotOutputBase_ == "runbased") {
0337 LuminosityBlockNumber_t firstLumi = beamSpotMap_.begin()->first;
0338 std::pair<edm::Timestamp, reco::BeamSpot> thepair(beamSpotMap_.begin()->second.first,
0339 weight(beamSpotMap_.begin(), beamSpotMap_.end()));
0340 beamSpotMap_.clear();
0341 beamSpotMap_[firstLumi] = thepair;
0342 } else {
0343 LogInfo("AlcaBeamSpotManager") << "Unrecognized BeamSpotOutputBase parameter: " << beamSpotOutputBase_ << endl;
0344 }
0345 }
0346
0347
0348 BeamSpot AlcaBeamSpotManager::weight(const bsMap_iterator &begin, const bsMap_iterator &end) {
0349 double x = 0, xError = 0;
0350 double y = 0, yError = 0;
0351 double z = 0, zError = 0;
0352 double sigmaZ = 0, sigmaZError = 0;
0353 double dxdz = 0, dxdzError = 0;
0354 double dydz = 0, dydzError = 0;
0355 double widthX = 0, widthXError = 0;
0356 double widthY = 0, widthYError = 0;
0357 LogInfo("AlcaBeamSpotManager") << "Weighted BeamSpot will span lumi " << begin->first << " to " << end->first << endl;
0358
0359 BeamSpot::BeamType type = BeamSpot::Unknown;
0360 for (bsMap_iterator it = begin; it != end; it++) {
0361 weight(x, xError, it->second.second.x0(), it->second.second.x0Error());
0362 weight(y, yError, it->second.second.y0(), it->second.second.y0Error());
0363 weight(z, zError, it->second.second.z0(), it->second.second.z0Error());
0364 weight(sigmaZ, sigmaZError, it->second.second.sigmaZ(), it->second.second.sigmaZ0Error());
0365 weight(dxdz, dxdzError, it->second.second.dxdz(), it->second.second.dxdzError());
0366 weight(dydz, dydzError, it->second.second.dydz(), it->second.second.dydzError());
0367 weight(widthX, widthXError, it->second.second.BeamWidthX(), it->second.second.BeamWidthXError());
0368 weight(widthY, widthYError, it->second.second.BeamWidthY(), it->second.second.BeamWidthYError());
0369 if (it->second.second.type() == BeamSpot::Tracker) {
0370 type = BeamSpot::Tracker;
0371 }
0372 }
0373 BeamSpot::Point bsPosition(x, y, z);
0374 BeamSpot::CovarianceMatrix error;
0375 error(0, 0) = xError * xError;
0376 error(1, 1) = yError * yError;
0377 error(2, 2) = zError * zError;
0378 error(3, 3) = sigmaZError * sigmaZError;
0379 error(4, 4) = dxdzError * dxdzError;
0380 error(5, 5) = dydzError * dydzError;
0381 error(6, 6) = widthXError * widthXError;
0382 BeamSpot weightedBeamSpot(bsPosition, sigmaZ, dxdz, dydz, widthX, error, type);
0383 weightedBeamSpot.setBeamWidthY(widthY);
0384 LogInfo("AlcaBeamSpotManager") << "Weighted BeamSpot will be:" << '\n' << weightedBeamSpot << endl;
0385 return weightedBeamSpot;
0386 }
0387
0388
0389 void AlcaBeamSpotManager::weight(double &mean, double &meanError, const double &val, const double &valError) {
0390 double tmpError = 0;
0391 if (meanError < 1e-8) {
0392 tmpError = 1 / (valError * valError);
0393 mean = val * tmpError;
0394 } else {
0395 tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
0396 mean = mean / (meanError * meanError) + val / (valError * valError);
0397 }
0398 mean = mean / tmpError;
0399 meanError = sqrt(1 / tmpError);
0400 }
0401
0402
0403 pair<float, float> AlcaBeamSpotManager::delta(const float &x,
0404 const float &xError,
0405 const float &nextX,
0406 const float &nextXError) {
0407 return pair<float, float>(x - nextX, sqrt(pow(xError, 2) + pow(nextXError, 2)));
0408 }
0409
0410
0411 float AlcaBeamSpotManager::deltaSig(const float &num, const float &den) {
0412 if (den != 0) {
0413 return fabs(num / den);
0414 } else {
0415 return float(LONG_MAX);
0416 }
0417 }
