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"))
54 aBeamSpot = beamSpotHandle.
product();
58 << *aBeamSpot << std::endl;
64 <<
" BS is not valid!" << std::endl;
71 vector<bsMap_iterator> listToErase;
74 listToErase.push_back(it);
77 for(vector<bsMap_iterator>::iterator it=listToErase.begin(); it !=listToErase.end(); it++){
95 map<LuminosityBlockNumber_t,BeamSpot> tmpBeamSpotMap_;
97 bool endOfRun =
false;
99 bool foundShift =
false;
111 <<
"Iteration: " << iteration <<
" size: " <<
beamSpotMap_.size() <<
"\n"
112 <<
"Lumi: " << currentBS->first <<
"\n"
114 <<
"\n" << nextBS->first <<
"\n" << nextBS->second
118 << nextNextBS->first <<
"\n" << nextNextBS->second
129 <<
"Reached lumi " << currentBS->first
130 <<
" now close payload because end of data has been reached.";
141 if (countlumi == maxNlumis -1){
143 <<
"close payload because maximum lumi sections accumulated within run ";
153 <<
"Checking checking!" << endl;
155 pair<float,float> adelta1;
156 pair<float,float> adelta2;
157 pair<float,float> adelta;
158 pair<float,float> adelta1dxdz;
159 pair<float,float> adelta2dxdz;
160 pair<float,float> adelta1dydz;
161 pair<float,float> adelta2dydz;
162 pair<float,float> adelta1widthX;
163 pair<float,float> adelta2widthX;
164 pair<float,float> adelta1widthY;
165 pair<float,float> adelta2widthY;
166 pair<float,float> adelta1z0;
167 pair<float,float> adelta1sigmaZ;
170 float min_limit = 0.0025;
173 limit = currentBS->second.BeamWidthX()/2.;
174 if(limit < min_limit){
179 adelta1 =
delta(currentBS->second.x0(), currentBS->second.x0Error(), nextBS->second.x0(), nextBS->second.x0Error());
180 adelta2 = pair<float,float>(0.,1.e9);
181 if (nextNextBS->second.type() != -1){
182 adelta2 =
delta(nextBS->second.x0(), nextBS->second.x0Error(), nextNextBS->second.x0(), nextNextBS->second.x0Error());
184 bool deltaX = (
deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >=
limit)?
true:
false;
186 if( !deltaX && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >=
limit){
188 <<
" positive, " << (adelta1.first+adelta2.first) <<
" limit=" << limit << endl;
191 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 <<
" negative, " << adelta1.first/adelta2.first << endl;
199 adelta1dxdz =
delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());
200 adelta2dxdz = pair<float,float>(0.,1.e9);
201 adelta1dydz =
delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());
202 adelta2dydz = pair<float,float>(0.,1.e9);
203 adelta1widthX =
delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());
204 adelta2widthX = pair<float,float>(0.,1.e9);
205 adelta1widthY =
delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());
206 adelta2widthY = pair<float,float>(0.,1.e9);
207 adelta1z0 =
delta(currentBS->second.z0(), currentBS->second.z0Error(), nextBS->second.z0(), nextBS->second.z0Error());
208 adelta1sigmaZ =
delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());
211 adelta1 =
delta(currentBS->second.y0(), currentBS->second.y0Error(), nextBS->second.y0(), nextBS->second.y0Error());
212 adelta2 = pair<float,float>(0.,1.e9);
214 adelta2 =
delta(nextBS->second.y0(), nextBS->second.y0Error(), nextNextBS->second.y0(), nextNextBS->second.y0Error());
215 adelta2dxdz =
delta(nextBS->second.dxdz(), nextBS->second.dxdzError(), nextNextBS->second.dxdz(), nextNextBS->second.dxdzError());
216 adelta2dydz =
delta(nextBS->second.dydz(), nextBS->second.dydzError(), nextNextBS->second.dydz(), nextNextBS->second.dydzError());
217 adelta2widthX =
delta(nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError(), nextNextBS->second.BeamWidthX(), nextNextBS->second.BeamWidthXError());
218 adelta2widthY =
delta(nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError(), nextNextBS->second.BeamWidthY(), nextNextBS->second.BeamWidthYError());
221 bool deltaY = (
deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >=
limit)?
true:
false;
223 if( !deltaY && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >=
limit){
225 <<
" positive, " << (adelta1.first+adelta2.first) <<
" limit=" << limit << endl;
228 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 <<
" negative, " << adelta1.first/adelta2.first << endl;
235 limit = currentBS->second.sigmaZ()/2.;
236 bool deltaZ = (
deltaSig(adelta1z0.first,adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >=
limit)?
true:
false;
237 adelta =
delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());
238 bool deltasigmaZ = (
deltaSig(adelta.first,adelta.second) > 5.0)?
true:
false;
239 bool deltadxdz =
false;
240 bool deltadydz =
false;
241 bool deltawidthX =
false;
242 bool deltawidthY =
false;
246 adelta =
delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());
247 deltadxdz = (
deltaSig(adelta.first,adelta.second) > 5.0)?
true:
false;
248 if(deltadxdz && (adelta1dxdz.first*adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 && fabs(adelta1dxdz.first/adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first/adelta2dxdz.first) < 3){
252 adelta =
delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());
253 deltadydz = (
deltaSig(adelta.first,adelta.second) > 5.0)?
true:
false;
254 if(deltadydz && (adelta1dydz.first*adelta2dydz.first) < 0 && adelta2dydz.first != 0 && fabs(adelta1dydz.first/adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first/adelta2dydz.first) < 3){
258 adelta =
delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());
259 deltawidthX = (
deltaSig(adelta.first,adelta.second) > 5.0)?
true:
false;
260 if(deltawidthX && (adelta1widthX.first*adelta2widthX.first) < 0 && adelta2widthX.first != 0 && fabs(adelta1widthX.first/adelta2widthX.first) > 0.33 && fabs(adelta1widthX.first/adelta2widthX.first) < 3){
264 adelta =
delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());
265 deltawidthY = (
deltaSig(adelta.first,adelta.second) > 5.0)?
true:
false;
266 if(deltawidthY && (adelta1widthY.first*adelta2widthY.first) < 0 && adelta2widthY.first != 0 && fabs(adelta1widthY.first/adelta2widthY.first) > 0.33 && fabs(adelta1widthY.first/adelta2widthY.first) < 3){
271 if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY){
275 <<
"close payload because of movement in"
279 <<
", sigmaZ=" << deltasigmaZ
280 <<
", dxdz=" << deltadxdz
281 <<
", dydz=" << deltadydz
282 <<
", widthX=" << deltawidthX
283 <<
", widthY=" << deltawidthY
290 tmpBeamSpotMap_[firstBS->first] =
weight(firstBS,nextBS);
294 tmpBeamSpotMap_[nextBS->first] = nextBS->second;
297 else if(!foundShift && !endOfRun){
298 tmpBeamSpotMap_[firstBS->first] =
weight(firstBS,nextBS);
337 double sigmaZ,sigmaZError = 0;
338 double dxdz,dxdzError = 0;
339 double dydz,dydzError = 0;
340 double widthX,widthXError = 0;
341 double widthY,widthYError = 0;
343 <<
"Weighted BeamSpot will span lumi "
344 << begin->first <<
" to " << end->first
349 weight(x , xError , it->second.x0() , it->second.x0Error());
350 weight(y , yError , it->second.y0() , it->second.y0Error());
351 weight(z , zError , it->second.z0() , it->second.z0Error());
352 weight(sigmaZ, sigmaZError, it->second.sigmaZ() , it->second.sigmaZ0Error());
353 weight(dxdz , dxdzError , it->second.dxdz() , it->second.dxdzError());
354 weight(dydz , dydzError , it->second.dydz() , it->second.dydzError());
355 weight(widthX, widthXError, it->second.BeamWidthX(), it->second.BeamWidthXError());
356 weight(widthY, widthYError, it->second.BeamWidthY(), it->second.BeamWidthYError());
363 error(0,0) = xError*xError;
364 error(1,1) = yError*yError;
365 error(2,2) = zError*zError;
366 error(3,3) = sigmaZError*sigmaZError;
367 error(4,4) = dxdzError*dxdzError;
368 error(5,5) = dydzError*dydzError;
369 error(6,6) = widthXError*widthXError;
370 BeamSpot weightedBeamSpot(bsPosition,sigmaZ,dxdz,dydz,widthX,error,type);
373 <<
"Weighted BeamSpot will be:" <<
'\n'
376 return weightedBeamSpot;
382 if (meanError < 1
e-8){
383 tmpError = 1/(valError*valError);
387 tmpError = 1/(meanError*meanError) + 1/(valError*valError);
388 mean = mean/(meanError*meanError) + val/(valError*valError);
390 mean = mean/tmpError;
391 meanError =
sqrt(1/tmpError);
396 return pair<float,float>(x - nextX,
sqrt(
pow(xError,2) +
pow(nextXError,2)) );
402 return fabs(num/den);
405 return float(LONG_MAX);
math::Error< dimension >::type CovarianceMatrix
std::string beamSpotOutputBase_
edm::InputTag beamSpotTag_
edm::EDGetTokenT< reco::BeamSpot > beamSpotToken_
AlcaBeamSpotManager(void)
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
unsigned int LuminosityBlockNumber_t
std::string beamSpotModuleName_
LuminosityBlockNumber_t luminosityBlock() const
void setBeamWidthY(double v)
reco::BeamSpot weight(const bsMap_iterator &begin, const bsMap_iterator &end)
virtual ~AlcaBeamSpotManager(void)
void createWeightedPayloads(void)
void readLumi(const edm::LuminosityBlock &)
T const * product() const
std::string beamSpotLabel_
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)
std::map< edm::LuminosityBlockNumber_t, reco::BeamSpot > beamSpotMap_
float deltaSig(const float &num, const float &den)