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"))
55 aBeamSpot = beamSpotHandle.
product();
59 << *aBeamSpot << std::endl;
65 <<
" BS is not valid!" << std::endl;
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);
310 if (!docreate) ++countlumi;
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;
345 <<
"Weighted BeamSpot will span lumi " 346 << begin->first <<
" to " << end->first
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());
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);
375 <<
"Weighted BeamSpot will be:" <<
'\n' 378 return weightedBeamSpot;
384 if (meanError < 1
e-8){
385 tmpError = 1/(valError*valError);
389 tmpError = 1/(meanError*meanError) + 1/(valError*valError);
390 mean = mean/(meanError*meanError) + val/(valError*valError);
392 mean = mean/tmpError;
393 meanError =
sqrt(1/tmpError);
398 return pair<float,float>(x - nextX,
sqrt(
pow(xError,2) +
pow(nextXError,2)) );
404 return fabs(num/den);
407 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)