23 : beamSpotOutputBase_(iConfig.getParameter<
ParameterSet>(
"AlcaBeamSpotHarvesterParameters")
24 .getUntrackedParameter<
std::
string>(
"BeamSpotOutputBase")),
25 beamSpotModuleName_(iConfig.getParameter<
ParameterSet>(
"AlcaBeamSpotHarvesterParameters")
26 .getUntrackedParameter<
std::
string>(
"BeamSpotModuleName")),
27 beamSpotLabel_(iConfig.getParameter<
ParameterSet>(
"AlcaBeamSpotHarvesterParameters")
28 .getUntrackedParameter<
std::
string>(
"BeamSpotLabel")),
29 sigmaZCut_(iConfig.getParameter<
ParameterSet>(
"AlcaBeamSpotHarvesterParameters")
30 .getUntrackedParameter<double>(
"SigmaZCut")) {
48 std::pair<edm::Timestamp, reco::BeamSpot> time_bs(iLumi.
beginTime(), *beamSpotHandle);
51 aBeamSpot = beamSpotHandle.
product();
53 LogInfo(
"AlcaBeamSpotManager") << *aBeamSpot << std::endl;
56 LogInfo(
"AlcaBeamSpotManager") <<
" BS is not valid!" << std::endl;
62 vector<bsMap_iterator> listToErase;
65 listToErase.push_back(it);
68 for (vector<bsMap_iterator>::iterator it = listToErase.begin(); it != listToErase.end(); it++) {
90 map<LuminosityBlockNumber_t, std::pair<edm::Timestamp, BeamSpot>> tmpBeamSpotMap_;
92 bool endOfRun =
false;
94 bool foundShift =
false;
105 <<
"Lumi: " << currentBS->first <<
"\n"
106 << currentBS->second.second <<
"\n"
107 << nextBS->first <<
"\n"
108 << nextBS->second.second << endl;
110 LogInfo(
"AlcaBeamSpotManager") << nextNextBS->first <<
"\n" << nextNextBS->second.second << endl;
118 LogInfo(
"AlcaBeamSpotManager") <<
"Reached lumi " << currentBS->first
119 <<
" now close payload because end of data has been reached.";
130 if (countlumi == maxNlumis - 1) {
131 LogInfo(
"AlcaBeamSpotManager") <<
"close payload because maximum lumi "
132 "sections accumulated within run ";
141 LogInfo(
"AlcaBeamSpotManager") <<
"Checking checking!" << endl;
143 pair<float, float> adelta1;
144 pair<float, float> adelta2;
145 pair<float, float> adelta;
146 pair<float, float> adelta1dxdz;
147 pair<float, float> adelta2dxdz;
148 pair<float, float> adelta1dydz;
149 pair<float, float> adelta2dydz;
150 pair<float, float> adelta1widthX;
151 pair<float, float> adelta2widthX;
152 pair<float, float> adelta1widthY;
153 pair<float, float> adelta2widthY;
154 pair<float, float> adelta1z0;
155 pair<float, float> adelta1sigmaZ;
158 float min_limit = 0.0025;
161 limit = currentBS->second.second.BeamWidthX() / 2.;
162 if (
limit < min_limit) {
166 currentBSObj = currentBS->second.second;
167 nextBSObj = nextBS->second.second;
171 adelta2 = pair<float, float>(0., 1.e9);
172 if (nextNextBS->second.second.type() != -1) {
173 adelta2 =
delta(nextBS->second.second.x0(),
175 nextNextBS->second.second.x0(),
176 nextNextBS->second.second.x0Error());
178 bool deltaX = (
deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >=
limit) ?
true :
false;
180 if (!deltaX && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >=
limit) {
182 <<
" positive, " << (adelta1.first + adelta2.first) <<
" limit=" <<
limit << endl;
184 }
else if (deltaX && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
185 fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
186 LogInfo(
"AlcaBeamSpotManager") <<
" negative, " << adelta1.first / adelta2.first << endl;
193 adelta2dxdz = pair<float, float>(0., 1.e9);
195 adelta2dydz = pair<float, float>(0., 1.e9);
200 adelta2widthX = pair<float, float>(0., 1.e9);
205 adelta2widthY = pair<float, float>(0., 1.e9);
212 adelta2 = pair<float, float>(0., 1.e9);
215 nextBSObj.
y0(), nextBSObj.
y0Error(), nextNextBS->second.second.y0(), nextNextBS->second.second.y0Error());
218 nextNextBS->second.second.dxdz(),
219 nextNextBS->second.second.dxdzError());
222 nextNextBS->second.second.dydz(),
223 nextNextBS->second.second.dydzError());
226 nextNextBS->second.second.BeamWidthX(),
227 nextNextBS->second.second.BeamWidthXError());
230 nextNextBS->second.second.BeamWidthY(),
231 nextNextBS->second.second.BeamWidthYError());
233 bool deltaY = (
deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >=
limit) ?
true :
false;
235 if (!deltaY && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >=
limit) {
237 <<
" positive, " << (adelta1.first + adelta2.first) <<
" limit=" <<
limit << endl;
239 }
else if (deltaY && adelta1.first * adelta2.first < 0 && adelta2.first != 0 &&
240 fabs(adelta1.first / adelta2.first) > 0.33 && fabs(adelta1.first / adelta2.first) < 3) {
241 LogInfo(
"AlcaBeamSpotManager") <<
" negative, " << adelta1.first / adelta2.first << endl;
248 (
deltaSig(adelta1z0.first, adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >=
limit) ?
true :
false;
251 bool deltasigmaZ = (
deltaSig(adelta.first, adelta.second) > 5.0) ?
true :
false;
252 bool deltadxdz =
false;
253 bool deltadydz =
false;
254 bool deltawidthX =
false;
255 bool deltawidthY =
false;
259 deltadxdz = (
deltaSig(adelta.first, adelta.second) > 5.0) ?
true :
false;
260 if (deltadxdz && (adelta1dxdz.first * adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 &&
261 fabs(adelta1dxdz.first / adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first / adelta2dxdz.first) < 3) {
266 deltadydz = (
deltaSig(adelta.first, adelta.second) > 5.0) ?
true :
false;
267 if (deltadydz && (adelta1dydz.first * adelta2dydz.first) < 0 && adelta2dydz.first != 0 &&
268 fabs(adelta1dydz.first / adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first / adelta2dydz.first) < 3) {
276 deltawidthX = (
deltaSig(adelta.first, adelta.second) > 5.0) ?
true :
false;
277 if (deltawidthX && (adelta1widthX.first * adelta2widthX.first) < 0 && adelta2widthX.first != 0 &&
278 fabs(adelta1widthX.first / adelta2widthX.first) > 0.33 &&
279 fabs(adelta1widthX.first / adelta2widthX.first) < 3) {
287 deltawidthY = (
deltaSig(adelta.first, adelta.second) > 5.0) ?
true :
false;
288 if (deltawidthY && (adelta1widthY.first * adelta2widthY.first) < 0 && adelta2widthY.first != 0 &&
289 fabs(adelta1widthY.first / adelta2widthY.first) > 0.33 &&
290 fabs(adelta1widthY.first / adelta2widthY.first) < 3) {
294 if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY) {
297 LogInfo(
"AlcaBeamSpotManager") <<
"close payload because of movement in"
298 <<
" X=" << deltaX <<
", Y=" << deltaY <<
", Z=" << deltaZ
299 <<
", sigmaZ=" << deltasigmaZ <<
", dxdz=" << deltadxdz
300 <<
", dydz=" << deltadydz <<
", widthX=" << deltawidthX
301 <<
", widthY=" << deltawidthY << endl;
305 std::pair<edm::Timestamp, reco::BeamSpot> thepair;
307 thepair = std::make_pair(currentBS->second.first,
weight(firstBS, nextBS));
308 tmpBeamSpotMap_[firstBS->first] = thepair;
312 thepair = std::make_pair(nextBS->second.first, nextBSObj);
313 tmpBeamSpotMap_[nextBS->first] = thepair;
315 }
else if (!foundShift && !endOfRun) {
316 thepair = std::make_pair(currentBS->second.first,
weight(firstBS, nextBS));
317 tmpBeamSpotMap_[firstBS->first] = thepair;
319 thepair = std::make_pair(currentBS->second.first,
weight(firstBS,
beamSpotMap_.end()));
320 tmpBeamSpotMap_[firstBS->first] = thepair;
339 std::pair<edm::Timestamp, reco::BeamSpot> thepair(
beamSpotMap_.begin()->second.first,
350 double x, xError = 0;
351 double y, yError = 0;
352 double z, zError = 0;
353 double sigmaZ, sigmaZError = 0;
354 double dxdz, dxdzError = 0;
355 double dydz, dydzError = 0;
356 double widthX, widthXError = 0;
357 double widthY, widthYError = 0;
358 LogInfo(
"AlcaBeamSpotManager") <<
"Weighted BeamSpot will span lumi " << begin->first <<
" to " <<
end->first << endl;
362 weight(
x, xError, it->second.second.x0(), it->second.second.x0Error());
363 weight(
y, yError, it->second.second.y0(), it->second.second.y0Error());
364 weight(
z, zError, it->second.second.z0(), it->second.second.z0Error());
365 weight(
sigmaZ, sigmaZError, it->second.second.sigmaZ(), it->second.second.sigmaZ0Error());
366 weight(
dxdz, dxdzError, it->second.second.dxdz(), it->second.second.dxdzError());
367 weight(
dydz, dydzError, it->second.second.dydz(), it->second.second.dydzError());
368 weight(
widthX, widthXError, it->second.second.BeamWidthX(), it->second.second.BeamWidthXError());
369 weight(
widthY, widthYError, it->second.second.BeamWidthY(), it->second.second.BeamWidthYError());
376 error(0, 0) = xError * xError;
377 error(1, 1) = yError * yError;
378 error(2, 2) = zError * zError;
379 error(3, 3) = sigmaZError * sigmaZError;
380 error(4, 4) = dxdzError * dxdzError;
381 error(5, 5) = dydzError * dydzError;
382 error(6, 6) = widthXError * widthXError;
385 LogInfo(
"AlcaBeamSpotManager") <<
"Weighted BeamSpot will be:" <<
'\n' << weightedBeamSpot << endl;
386 return weightedBeamSpot;
392 if (meanError < 1
e-8) {
393 tmpError = 1 / (valError * valError);
396 tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
397 mean =
mean / (meanError * meanError) +
val / (valError * valError);
400 meanError =
sqrt(1 / tmpError);
407 const float &nextXError) {
408 return pair<float, float>(
x - nextX,
sqrt(
pow(xError, 2) +
pow(nextXError, 2)));
414 return fabs(
num / den);
416 return float(LONG_MAX);