CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
AlcaBeamSpotManager.cc
Go to the documentation of this file.
1 
10 #include <climits>
11 #include <cmath>
12 #include <vector>
13 
14 using namespace edm;
15 using namespace reco;
16 using namespace std;
17 
18 //--------------------------------------------------------------------------------------------------
20 
21 //--------------------------------------------------------------------------------------------------
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")) {
33  LogInfo("AlcaBeamSpotManager") << "Output base: " << beamSpotOutputBase_ << std::endl;
34  reset();
35 }
36 
37 //--------------------------------------------------------------------------------------------------
39 
40 //--------------------------------------------------------------------------------------------------
42 //--------------------------------------------------------------------------------------------------
44  Handle<BeamSpot> beamSpotHandle;
45  iLumi.getByToken(beamSpotToken_, beamSpotHandle);
46 
47  if (beamSpotHandle.isValid()) { // check the product
48  std::pair<edm::Timestamp, reco::BeamSpot> time_bs(iLumi.beginTime(), *beamSpotHandle);
49  beamSpotMap_[iLumi.luminosityBlock()] = time_bs;
50  const BeamSpot *aBeamSpot = &beamSpotMap_[iLumi.luminosityBlock()].second;
51  aBeamSpot = beamSpotHandle.product();
52  LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
53  LogInfo("AlcaBeamSpotManager") << *aBeamSpot << std::endl;
54  } else {
55  LogInfo("AlcaBeamSpotManager") << "Lumi: " << iLumi.luminosityBlock() << std::endl;
56  LogInfo("AlcaBeamSpotManager") << " BS is not valid!" << std::endl;
57  }
58 }
59 
60 //--------------------------------------------------------------------------------------------------
62  vector<bsMap_iterator> listToErase;
63  for (bsMap_iterator it = beamSpotMap_.begin(); it != beamSpotMap_.end(); it++) {
64  if (it->second.second.type() != BeamSpot::Tracker || it->second.second.sigmaZ() < sigmaZCut_) {
65  listToErase.push_back(it);
66  }
67  }
68  for (vector<bsMap_iterator>::iterator it = listToErase.begin(); it != listToErase.end(); it++) {
69  beamSpotMap_.erase(*it);
70  }
71  if (beamSpotMap_.size() <= 1) {
72  return;
73  }
74  // Return only if lumibased since the collapsing alghorithm requires the next
75  // and next to next lumi sections
76  else if (beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased") {
77  return;
78  }
79  if (beamSpotOutputBase_ == "lumibased") {
80  // bsMap_iterator referenceBS = beamSpotMap_.begin();
81  bsMap_iterator firstBS = beamSpotMap_.begin();
82  // bsMap_iterator lastBS = beamSpotMap_.begin();
83  bsMap_iterator currentBS = beamSpotMap_.begin();
84  bsMap_iterator nextBS = ++beamSpotMap_.begin();
85  bsMap_iterator nextNextBS = ++(++(beamSpotMap_.begin()));
86 
87  reco::BeamSpot currentBSObj;
88  reco::BeamSpot nextBSObj;
89 
90  map<LuminosityBlockNumber_t, std::pair<edm::Timestamp, BeamSpot>> tmpBeamSpotMap_;
91  bool docreate = true;
92  bool endOfRun = false; // Added
93  bool docheck = true; // Added
94  bool foundShift = false;
95  long countlumi = 0; // Added
96  string tmprun = ""; // Added
97  long maxNlumis = 20; // Added
98  // if weighted:
99  // maxNlumis = 999999999
100 
101  unsigned int iteration = 0;
102  // while(nextNextBS!=beamSpotMap_.end()){
103  while (nextBS != beamSpotMap_.end()) {
104  LogInfo("AlcaBeamSpotManager") << "Iteration: " << iteration << " size: " << beamSpotMap_.size() << "\n"
105  << "Lumi: " << currentBS->first << "\n"
106  << currentBS->second.second << "\n"
107  << nextBS->first << "\n"
108  << nextBS->second.second << endl;
109  if (nextNextBS != beamSpotMap_.end())
110  LogInfo("AlcaBeamSpotManager") << nextNextBS->first << "\n" << nextNextBS->second.second << endl;
111 
112  if (docreate) {
113  firstBS = currentBS;
114  docreate = false; // Added
115  }
116  // if(iteration >= beamSpotMap_.size()-3){
117  if (iteration >= beamSpotMap_.size() - 2) {
118  LogInfo("AlcaBeamSpotManager") << "Reached lumi " << currentBS->first
119  << " now close payload because end of data has been reached.";
120  docreate = true;
121  endOfRun = true;
122  }
123  // check we run over the same run
124  // if (ibeam->first.Run() != inextbeam->first.Run()){
125  // LogInfo("AlcaBeamSpotManager")
126  // << "close payload because end of run.";
127  // docreate = true;
128  // }
129  // check maximum lumi counts
130  if (countlumi == maxNlumis - 1) {
131  LogInfo("AlcaBeamSpotManager") << "close payload because maximum lumi "
132  "sections accumulated within run ";
133  docreate = true;
134  countlumi = 0;
135  }
136  // if weighted:
137  // docheck = False
138  // check offsets
139  if (docheck) {
140  foundShift = false;
141  LogInfo("AlcaBeamSpotManager") << "Checking checking!" << endl;
142  float limit = 0;
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;
156 
157  // define minimum limit
158  float min_limit = 0.0025;
159 
160  // limit for x and y
161  limit = currentBS->second.second.BeamWidthX() / 2.;
162  if (limit < min_limit) {
163  limit = min_limit;
164  }
165 
166  currentBSObj = currentBS->second.second;
167  nextBSObj = nextBS->second.second;
168 
169  // check movements in X
170  adelta1 = delta(currentBSObj.x0(), currentBSObj.x0Error(), nextBSObj.x0(), nextBSObj.x0Error());
171  adelta2 = pair<float, float>(0., 1.e9);
172  if (nextNextBS->second.second.type() != -1) {
173  adelta2 = delta(nextBS->second.second.x0(),
174  nextBSObj.x0Error(),
175  nextNextBS->second.second.x0(),
176  nextNextBS->second.second.x0Error());
177  }
178  bool deltaX = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
179  if (iteration < beamSpotMap_.size() - 2) {
180  if (!deltaX && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
181  LogInfo("AlcaBeamSpotManager")
182  << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
183  deltaX = true;
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;
187  deltaX = false;
188  }
189  }
190 
191  // calculating all deltas
192  adelta1dxdz = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
193  adelta2dxdz = pair<float, float>(0., 1.e9);
194  adelta1dydz = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
195  adelta2dydz = pair<float, float>(0., 1.e9);
196  adelta1widthX = delta(currentBSObj.BeamWidthX(),
197  currentBSObj.BeamWidthXError(),
198  nextBSObj.BeamWidthX(),
199  nextBSObj.BeamWidthXError());
200  adelta2widthX = pair<float, float>(0., 1.e9);
201  adelta1widthY = delta(currentBSObj.BeamWidthY(),
202  currentBSObj.BeamWidthYError(),
203  nextBSObj.BeamWidthY(),
204  nextBSObj.BeamWidthYError());
205  adelta2widthY = pair<float, float>(0., 1.e9);
206  adelta1z0 = delta(currentBSObj.z0(), currentBSObj.z0Error(), nextBSObj.z0(), nextBSObj.z0Error());
207  adelta1sigmaZ =
208  delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
209 
210  // check movements in Y
211  adelta1 = delta(currentBSObj.y0(), currentBSObj.y0Error(), nextBSObj.y0(), nextBSObj.y0Error());
212  adelta2 = pair<float, float>(0., 1.e9);
213  if (nextNextBS->second.second.type() != BeamSpot::Unknown) {
214  adelta2 = delta(
215  nextBSObj.y0(), nextBSObj.y0Error(), nextNextBS->second.second.y0(), nextNextBS->second.second.y0Error());
216  adelta2dxdz = delta(nextBSObj.dxdz(),
217  nextBSObj.dxdzError(),
218  nextNextBS->second.second.dxdz(),
219  nextNextBS->second.second.dxdzError());
220  adelta2dydz = delta(nextBSObj.dydz(),
221  nextBSObj.dydzError(),
222  nextNextBS->second.second.dydz(),
223  nextNextBS->second.second.dydzError());
224  adelta2widthX = delta(nextBSObj.BeamWidthX(),
225  nextBSObj.BeamWidthXError(),
226  nextNextBS->second.second.BeamWidthX(),
227  nextNextBS->second.second.BeamWidthXError());
228  adelta2widthY = delta(nextBSObj.BeamWidthY(),
229  nextBSObj.BeamWidthYError(),
230  nextNextBS->second.second.BeamWidthY(),
231  nextNextBS->second.second.BeamWidthYError());
232  }
233  bool deltaY = (deltaSig(adelta1.first, adelta1.second) > 3.5 && adelta1.first >= limit) ? true : false;
234  if (iteration < beamSpotMap_.size() - 2) {
235  if (!deltaY && adelta1.first * adelta2.first > 0. && fabs(adelta1.first + adelta2.first) >= limit) {
236  LogInfo("AlcaBeamSpotManager")
237  << " positive, " << (adelta1.first + adelta2.first) << " limit=" << limit << endl;
238  deltaY = true;
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;
242  deltaY = false;
243  }
244  }
245 
246  limit = currentBSObj.sigmaZ() / 2.;
247  bool deltaZ =
248  (deltaSig(adelta1z0.first, adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >= limit) ? true : false;
249  adelta =
250  delta(currentBSObj.sigmaZ(), currentBSObj.sigmaZ0Error(), nextBSObj.sigmaZ(), nextBSObj.sigmaZ0Error());
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;
256 
257  if (iteration < beamSpotMap_.size() - 2) {
258  adelta = delta(currentBSObj.dxdz(), currentBSObj.dxdzError(), nextBSObj.dxdz(), nextBSObj.dxdzError());
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) {
262  deltadxdz = false;
263  }
264 
265  adelta = delta(currentBSObj.dydz(), currentBSObj.dydzError(), nextBSObj.dydz(), nextBSObj.dydzError());
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) {
269  deltadydz = false;
270  }
271 
272  adelta = delta(currentBSObj.BeamWidthX(),
273  currentBSObj.BeamWidthXError(),
274  nextBSObj.BeamWidthX(),
275  nextBSObj.BeamWidthXError());
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) {
280  deltawidthX = false;
281  }
282 
283  adelta = delta(currentBSObj.BeamWidthY(),
284  currentBSObj.BeamWidthYError(),
285  nextBSObj.BeamWidthY(),
286  nextBSObj.BeamWidthYError());
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) {
291  deltawidthY = false;
292  }
293  }
294  if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY) {
295  docreate = true;
296  foundShift = true;
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;
302  }
303  }
304  if (docreate) {
305  std::pair<edm::Timestamp, reco::BeamSpot> thepair;
306  if (foundShift) {
307  thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
308  tmpBeamSpotMap_[firstBS->first] = thepair;
309  if (endOfRun) {
310  // if we're here, then we need to found a shift in the last LS
311  // We already created a new IOV, now create one just for the last LS
312  thepair = std::make_pair(nextBS->second.first, nextBSObj);
313  tmpBeamSpotMap_[nextBS->first] = thepair;
314  }
315  } else if (!foundShift && !endOfRun) { // maxLS reached
316  thepair = std::make_pair(currentBS->second.first, weight(firstBS, nextBS));
317  tmpBeamSpotMap_[firstBS->first] = thepair;
318  } else { // end of run with no shift detectred in last LS
319  thepair = std::make_pair(currentBS->second.first, weight(firstBS, beamSpotMap_.end()));
320  tmpBeamSpotMap_[firstBS->first] = thepair;
321  }
322  firstBS = nextBS;
323  countlumi = 0;
324  }
325  // tmprun = currentBS->second.Run
326  // increase the counter by one only if the IOV hasn't been closed
327  if (!docreate)
328  ++countlumi;
329 
330  currentBS = nextBS;
331  nextBS = nextNextBS;
332  nextNextBS++;
333  ++iteration;
334  }
335  beamSpotMap_.clear();
336  beamSpotMap_ = tmpBeamSpotMap_;
337  } else if (beamSpotOutputBase_ == "runbased") {
338  LuminosityBlockNumber_t firstLumi = beamSpotMap_.begin()->first;
339  std::pair<edm::Timestamp, reco::BeamSpot> thepair(beamSpotMap_.begin()->second.first,
340  weight(beamSpotMap_.begin(), beamSpotMap_.end()));
341  beamSpotMap_.clear();
342  beamSpotMap_[firstLumi] = thepair;
343  } else {
344  LogInfo("AlcaBeamSpotManager") << "Unrecognized BeamSpotOutputBase parameter: " << beamSpotOutputBase_ << endl;
345  }
346 }
347 
348 //--------------------------------------------------------------------------------------------------
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;
359 
361  for (bsMap_iterator it = begin; it != end; it++) {
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());
370  if (it->second.second.type() == BeamSpot::Tracker) {
371  type = BeamSpot::Tracker;
372  }
373  }
374  BeamSpot::Point bsPosition(x, y, z);
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;
383  BeamSpot weightedBeamSpot(bsPosition, sigmaZ, dxdz, dydz, widthX, error, type);
384  weightedBeamSpot.setBeamWidthY(widthY);
385  LogInfo("AlcaBeamSpotManager") << "Weighted BeamSpot will be:" << '\n' << weightedBeamSpot << endl;
386  return weightedBeamSpot;
387 }
388 
389 //--------------------------------------------------------------------------------------------------
390 void AlcaBeamSpotManager::weight(double &mean, double &meanError, const double &val, const double &valError) {
391  double tmpError = 0;
392  if (meanError < 1e-8) {
393  tmpError = 1 / (valError * valError);
394  mean = val * tmpError;
395  } else {
396  tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
397  mean = mean / (meanError * meanError) + val / (valError * valError);
398  }
399  mean = mean / tmpError;
400  meanError = sqrt(1 / tmpError);
401 }
402 
403 //--------------------------------------------------------------------------------------------------
404 pair<float, float> AlcaBeamSpotManager::delta(const float &x,
405  const float &xError,
406  const float &nextX,
407  const float &nextXError) {
408  return pair<float, float>(x - nextX, sqrt(pow(xError, 2) + pow(nextXError, 2)));
409 }
410 
411 //--------------------------------------------------------------------------------------------------
412 float AlcaBeamSpotManager::deltaSig(const float &num, const float &den) {
413  if (den != 0) {
414  return fabs(num / den);
415  } else {
416  return float(LONG_MAX);
417  }
418 }
math::Error< dimension >::type CovarianceMatrix
Definition: BeamSpot.h:29
double z0() const
z coordinate
Definition: BeamSpot.h:65
BeamType
beam spot flags
Definition: BeamSpot.h:24
double sigmaZ0Error() const
error on sigma z
Definition: BeamSpot.h:92
std::string beamSpotOutputBase_
std::map< edm::LuminosityBlockNumber_t, std::pair< edm::Timestamp, reco::BeamSpot > > beamSpotMap_
edm::InputTag beamSpotTag_
edm::EDGetTokenT< reco::BeamSpot > beamSpotToken_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
double dydzError() const
error on dydz
Definition: BeamSpot.h:96
std::map< edm::LuminosityBlockNumber_t, std::pair< edm::Timestamp, reco::BeamSpot > >::iterator bsMap_iterator
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:27
Timestamp const & beginTime() const
unsigned int LuminosityBlockNumber_t
std::string beamSpotModuleName_
LuminosityBlockNumber_t luminosityBlock() const
double dydz() const
dydz slope
Definition: BeamSpot.h:80
tuple iteration
Definition: align_cfg.py:5
void setBeamWidthY(double v)
Definition: BeamSpot.h:105
double dxdzError() const
error on dxdz
Definition: BeamSpot.h:94
T sqrt(T t)
Definition: SSEVec.h:19
reco::BeamSpot weight(const bsMap_iterator &begin, const bsMap_iterator &end)
uint16_t const *__restrict__ x
Definition: gpuClustering.h:39
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:82
virtual ~AlcaBeamSpotManager(void)
double BeamWidthYError() const
error on beam width Y, assume error in X = Y
Definition: BeamSpot.h:101
bool isValid() const
Definition: HandleBase.h:70
double BeamWidthXError() const
error on beam width X, assume error in X = Y
Definition: BeamSpot.h:99
double z0Error() const
error on z
Definition: BeamSpot.h:90
double dxdz() const
dxdz slope
Definition: BeamSpot.h:78
double x0Error() const
error on x
Definition: BeamSpot.h:86
double y0Error() const
error on y
Definition: BeamSpot.h:88
tuple tmprun
Definition: ntuplemaker.py:244
Log< level::Info, false > LogInfo
void readLumi(const edm::LuminosityBlock &)
T const * product() const
Definition: Handle.h:70
double sigmaZ() const
sigma z
Definition: BeamSpot.h:76
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:84
string end
Definition: dataset.py:937
double y0() const
y coordinate
Definition: BeamSpot.h:63
std::pair< float, float > delta(const float &x, const float &xError, const float &nextX, const float &nextXError)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
float deltaSig(const float &num, const float &den)
double x0() const
x coordinate
Definition: BeamSpot.h:61