CMS 3D CMS Logo

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