File indexing completed on 2024-10-12 04:20:21
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 = &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
0075
0076 else if (beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased") {
0077 return;
0078 }
0079 if (beamSpotOutputBase_ == "lumibased") {
0080
0081 bsMap_iterator firstBS = beamSpotMap_.begin();
0082
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;
0093 bool docheck = true;
0094 bool foundShift = false;
0095 long countlumi = 0;
0096 string tmprun = "";
0097 long maxNlumis = 20;
0098
0099
0100
0101 unsigned int iteration = 0;
0102
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;
0115 }
0116
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
0124
0125
0126
0127
0128
0129
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
0137
0138
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
0158 float min_limit = 0.0025;
0159
0160
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
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
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
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
0311
0312 thepair = std::make_pair(nextBS->second.first, nextBSObj);
0313 tmpBeamSpotMap_[nextBS->first] = thepair;
0314 }
0315 } else if (!foundShift && !endOfRun) {
0316 thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
0317 tmpBeamSpotMap_[firstBS->first] = thepair;
0318 } else {
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
0326
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 = 0, xError = 0;
0351 double y = 0, yError = 0;
0352 double z = 0, zError = 0;
0353 double sigmaZ = 0, sigmaZError = 0;
0354 double dxdz = 0, dxdzError = 0;
0355 double dydz = 0, dydzError = 0;
0356 double widthX = 0, widthXError = 0;
0357 double widthY = 0, 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 }
