72 vector<bsMap_iterator> listToErase;
75 listToErase.push_back(it);
78 for(vector<bsMap_iterator>::iterator it=listToErase.begin(); it !=listToErase.end(); it++){
96 map<LuminosityBlockNumber_t,BeamSpot> tmpBeamSpotMap_;
98 bool endOfRun =
false;
100 bool foundShift =
false;
112 <<
"Iteration: " << iteration <<
" size: " <<
beamSpotMap_.size() <<
"\n"
113 <<
"Lumi: " << currentBS->first <<
"\n"
115 <<
"\n" << nextBS->first <<
"\n" << nextBS->second
119 << nextNextBS->first <<
"\n" << nextNextBS->second
130 <<
"Reached lumi " << currentBS->first
131 <<
" now close payload because end of data has been reached.";
142 if (countlumi == maxNlumis -1){
144 <<
"close payload because maximum lumi sections accumulated within run ";
154 <<
"Checking checking!" << endl;
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;
171 float min_limit = 0.0025;
174 limit = currentBS->second.BeamWidthX()/2.;
175 if(limit < min_limit){
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());
185 bool deltaX = (
deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >=
limit)?
true:
false;
187 if( !deltaX && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >=
limit){
189 <<
" positive, " << (adelta1.first+adelta2.first) <<
" limit=" << limit << endl;
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){
194 <<
" negative, " << adelta1.first/adelta2.first << endl;
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());
212 adelta1 =
delta(currentBS->second.y0(), currentBS->second.y0Error(), nextBS->second.y0(), nextBS->second.y0Error());
213 adelta2 = pair<float,float>(0.,1.e9);
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());
222 bool deltaY = (
deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >=
limit)?
true:
false;
224 if( !deltaY && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >=
limit){
226 <<
" positive, " << (adelta1.first+adelta2.first) <<
" limit=" << limit << endl;
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){
231 <<
" negative, " << adelta1.first/adelta2.first << endl;
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;
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){
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){
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){
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){
272 if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY){
276 <<
"close payload because of movement in"
280 <<
", sigmaZ=" << deltasigmaZ
281 <<
", dxdz=" << deltadxdz
282 <<
", dydz=" << deltadydz
283 <<
", widthX=" << deltawidthX
284 <<
", widthY=" << deltawidthY
291 tmpBeamSpotMap_[firstBS->first] =
weight(firstBS,nextBS);
295 tmpBeamSpotMap_[nextBS->first] = nextBS->second;
298 else if(!foundShift && !endOfRun){
299 tmpBeamSpotMap_[firstBS->first] =
weight(firstBS,nextBS);
std::string beamSpotOutputBase_
std::map< edm::LuminosityBlockNumber_t, reco::BeamSpot >::iterator bsMap_iterator
unsigned int LuminosityBlockNumber_t
reco::BeamSpot weight(const bsMap_iterator &begin, const bsMap_iterator &end)
std::pair< float, float > delta(const float &x, const float &xError, const float &nextX, const float &nextXError)
std::map< edm::LuminosityBlockNumber_t, reco::BeamSpot > beamSpotMap_
float deltaSig(const float &num, const float &den)