00001
00009 #include "Calibration/TkAlCaRecoProducers/interface/AlcaBeamSpotManager.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 #include "FWCore/Framework/interface/ESHandle.h"
00013 #include <vector>
00014 #include <math.h>
00015 #include <limits.h>
00016
00017 using namespace edm;
00018 using namespace reco;
00019 using namespace std;
00020
00021
00022 AlcaBeamSpotManager::AlcaBeamSpotManager(void){
00023 }
00024
00025
00026 AlcaBeamSpotManager::AlcaBeamSpotManager(const ParameterSet& iConfig) :
00027 beamSpotOutputBase_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotOutputBase")),
00028 beamSpotModuleName_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotModuleName")),
00029 beamSpotLabel_ (iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotLabel"))
00030 {
00031 LogInfo("AlcaBeamSpotManager")
00032 << "Output base: " << beamSpotOutputBase_
00033 << std::endl;
00034 reset();
00035 }
00036
00037
00038 AlcaBeamSpotManager::~AlcaBeamSpotManager(void){
00039 }
00040
00041
00042 void AlcaBeamSpotManager::reset(void){
00043 beamSpotMap_.clear();
00044 }
00045
00046 void AlcaBeamSpotManager::readLumi(const LuminosityBlock& iLumi){
00047
00048 Handle<BeamSpot> beamSpotHandle;
00049 iLumi.getByLabel(beamSpotModuleName_,beamSpotLabel_, beamSpotHandle);
00050
00051
00052 if(beamSpotHandle.isValid()) {
00053 beamSpotMap_[iLumi.luminosityBlock()] = *beamSpotHandle;
00054 const BeamSpot* aBeamSpot = &beamSpotMap_[iLumi.luminosityBlock()];
00055 aBeamSpot = beamSpotHandle.product();
00056 LogInfo("AlcaBeamSpotManager")
00057 << "Lumi: " << iLumi.luminosityBlock() << std::endl;
00058 LogInfo("AlcaBeamSpotManager")
00059 << *aBeamSpot << std::endl;
00060 }
00061 else {
00062 LogInfo("AlcaBeamSpotManager")
00063 << "Lumi: " << iLumi.luminosityBlock() << std::endl;
00064 LogInfo("AlcaBeamSpotManager")
00065 << " BS is not valid!" << std::endl;
00066 }
00067
00068 }
00069
00070
00071 void AlcaBeamSpotManager::createWeightedPayloads(void){
00072 vector<bsMap_iterator> listToErase;
00073 for(bsMap_iterator it=beamSpotMap_.begin(); it!=beamSpotMap_.end();it++){
00074 if(it->second.type() != BeamSpot::Tracker){
00075 listToErase.push_back(it);
00076 }
00077 }
00078 for(vector<bsMap_iterator>::iterator it=listToErase.begin(); it !=listToErase.end(); it++){
00079 beamSpotMap_.erase(*it);
00080 }
00081 if(beamSpotMap_.size() <= 1){
00082 return;
00083 }
00084
00085 else if(beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased"){
00086 return;
00087 }
00088 if(beamSpotOutputBase_ == "lumibased"){
00089
00090 bsMap_iterator firstBS = beamSpotMap_.begin();
00091
00092 bsMap_iterator currentBS = beamSpotMap_.begin();
00093 bsMap_iterator nextBS = ++beamSpotMap_.begin();
00094 bsMap_iterator nextNextBS = ++(++(beamSpotMap_.begin()));
00095
00096 map<LuminosityBlockNumber_t,BeamSpot> tmpBeamSpotMap_;
00097 bool docreate = true;
00098 bool endOfRun = false;
00099 bool docheck = true;
00100 bool foundShift = false;
00101 long countlumi = 0;
00102 string tmprun = "";
00103 long maxNlumis = 60;
00104
00105
00106
00107
00108 unsigned int iteration = 0;
00109
00110 while(nextBS!=beamSpotMap_.end()){
00111 LogInfo("AlcaBeamSpotManager")
00112 << "Iteration: " << iteration << " size: " << beamSpotMap_.size() << "\n"
00113 << "Lumi: " << currentBS->first << "\n"
00114 << currentBS->second
00115 << "\n" << nextBS->first << "\n" << nextBS->second
00116 << endl;
00117 if (nextNextBS!=beamSpotMap_.end())
00118 LogInfo("AlcaBeamSpotManager")
00119 << nextNextBS->first << "\n" << nextNextBS->second
00120 << endl;
00121
00122
00123 if(docreate){
00124 firstBS = currentBS;
00125 docreate = false;
00126 }
00127
00128 if(iteration >= beamSpotMap_.size()-2){
00129 LogInfo("AlcaBeamSpotManager")
00130 << "Reached lumi " << currentBS->first
00131 << " now close payload because end of data has been reached.";
00132 docreate = true;
00133 endOfRun = true;
00134 }
00135
00136
00137
00138
00139
00140
00141
00142 if (countlumi == maxNlumis -1){
00143 LogInfo("AlcaBeamSpotManager")
00144 << "close payload because maximum lumi sections accumulated within run ";
00145 docreate = true;
00146 countlumi = 0;
00147 }
00148
00149
00150
00151 if(docheck){
00152 foundShift = false;
00153 LogInfo("AlcaBeamSpotManager")
00154 << "Checking checking!" << endl;
00155 float limit = 0;
00156 pair<float,float> adelta1;
00157 pair<float,float> adelta2;
00158 pair<float,float> adelta;
00159 pair<float,float> adelta1dxdz;
00160 pair<float,float> adelta2dxdz;
00161 pair<float,float> adelta1dydz;
00162 pair<float,float> adelta2dydz;
00163 pair<float,float> adelta1widthX;
00164 pair<float,float> adelta2widthX;
00165 pair<float,float> adelta1widthY;
00166 pair<float,float> adelta2widthY;
00167 pair<float,float> adelta1z0;
00168 pair<float,float> adelta1sigmaZ;
00169
00170
00171 float min_limit = 0.0025;
00172
00173
00174 limit = currentBS->second.BeamWidthX()/2.;
00175 if(limit < min_limit){
00176 limit = min_limit;
00177 }
00178
00179
00180 adelta1 = delta(currentBS->second.x0(), currentBS->second.x0Error(), nextBS->second.x0(), nextBS->second.x0Error());
00181 adelta2 = pair<float,float>(0.,1.e9);
00182 if (nextNextBS->second.type() != -1){
00183 adelta2 = delta(nextBS->second.x0(), nextBS->second.x0Error(), nextNextBS->second.x0(), nextNextBS->second.x0Error());
00184 }
00185 bool deltaX = (deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >= limit)?true:false;
00186 if(iteration < beamSpotMap_.size()-2){
00187 if( !deltaX && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >= limit){
00188 LogInfo("AlcaBeamSpotManager")
00189 << " positive, " << (adelta1.first+adelta2.first) << " limit=" << limit << endl;
00190 deltaX = true;
00191 }
00192 else if( deltaX && adelta1.first*adelta2.first < 0 && adelta2.first != 0 && fabs(adelta1.first/adelta2.first) > 0.33 && fabs(adelta1.first/adelta2.first) < 3){
00193 LogInfo("AlcaBeamSpotManager")
00194 << " negative, " << adelta1.first/adelta2.first << endl;
00195 deltaX = false;
00196 }
00197 }
00198
00199
00200 adelta1dxdz = delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());
00201 adelta2dxdz = pair<float,float>(0.,1.e9);
00202 adelta1dydz = delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());
00203 adelta2dydz = pair<float,float>(0.,1.e9);
00204 adelta1widthX = delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());
00205 adelta2widthX = pair<float,float>(0.,1.e9);
00206 adelta1widthY = delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());
00207 adelta2widthY = pair<float,float>(0.,1.e9);
00208 adelta1z0 = delta(currentBS->second.z0(), currentBS->second.z0Error(), nextBS->second.z0(), nextBS->second.z0Error());
00209 adelta1sigmaZ = delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());
00210
00211
00212 adelta1 = delta(currentBS->second.y0(), currentBS->second.y0Error(), nextBS->second.y0(), nextBS->second.y0Error());
00213 adelta2 = pair<float,float>(0.,1.e9);
00214 if( nextNextBS->second.type() != BeamSpot::Unknown){
00215 adelta2 = delta(nextBS->second.y0(), nextBS->second.y0Error(), nextNextBS->second.y0(), nextNextBS->second.y0Error());
00216 adelta2dxdz = delta(nextBS->second.dxdz(), nextBS->second.dxdzError(), nextNextBS->second.dxdz(), nextNextBS->second.dxdzError());
00217 adelta2dydz = delta(nextBS->second.dydz(), nextBS->second.dydzError(), nextNextBS->second.dydz(), nextNextBS->second.dydzError());
00218 adelta2widthX = delta(nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError(), nextNextBS->second.BeamWidthX(), nextNextBS->second.BeamWidthXError());
00219 adelta2widthY = delta(nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError(), nextNextBS->second.BeamWidthY(), nextNextBS->second.BeamWidthYError());
00220
00221 }
00222 bool deltaY = (deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >= limit)?true:false;
00223 if(iteration < beamSpotMap_.size()-2){
00224 if( !deltaY && adelta1.first*adelta2.first > 0. && fabs(adelta1.first+adelta2.first) >= limit){
00225 LogInfo("AlcaBeamSpotManager")
00226 << " positive, " << (adelta1.first+adelta2.first) << " limit=" << limit << endl;
00227 deltaY = true;
00228 }
00229 else if( deltaY && adelta1.first*adelta2.first < 0 && adelta2.first != 0 && fabs(adelta1.first/adelta2.first) > 0.33 && fabs(adelta1.first/adelta2.first) < 3){
00230 LogInfo("AlcaBeamSpotManager")
00231 << " negative, " << adelta1.first/adelta2.first << endl;
00232 deltaY = false;
00233 }
00234 }
00235
00236 limit = currentBS->second.sigmaZ()/2.;
00237 bool deltaZ = (deltaSig(adelta1z0.first,adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >= limit)?true:false;
00238 adelta = delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());
00239 bool deltasigmaZ = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;
00240 bool deltadxdz = false;
00241 bool deltadydz = false;
00242 bool deltawidthX = false;
00243 bool deltawidthY = false;
00244
00245 if(iteration < beamSpotMap_.size()-2){
00246
00247 adelta = delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());
00248 deltadxdz = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;
00249 if(deltadxdz && (adelta1dxdz.first*adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 && fabs(adelta1dxdz.first/adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first/adelta2dxdz.first) < 3){
00250 deltadxdz = false;
00251 }
00252
00253 adelta = delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());
00254 deltadydz = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;
00255 if(deltadydz && (adelta1dydz.first*adelta2dydz.first) < 0 && adelta2dydz.first != 0 && fabs(adelta1dydz.first/adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first/adelta2dydz.first) < 3){
00256 deltadydz = false;
00257 }
00258
00259 adelta = delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());
00260 deltawidthX = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;
00261 if(deltawidthX && (adelta1widthX.first*adelta2widthX.first) < 0 && adelta2widthX.first != 0 && fabs(adelta1widthX.first/adelta2widthX.first) > 0.33 && fabs(adelta1widthX.first/adelta2widthX.first) < 3){
00262 deltawidthX = false;
00263 }
00264
00265 adelta = delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());
00266 deltawidthY = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;
00267 if(deltawidthY && (adelta1widthY.first*adelta2widthY.first) < 0 && adelta2widthY.first != 0 && fabs(adelta1widthY.first/adelta2widthY.first) > 0.33 && fabs(adelta1widthY.first/adelta2widthY.first) < 3){
00268 deltawidthY = false;
00269 }
00270
00271 }
00272 if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY){
00273 docreate = true;
00274 foundShift = true;
00275 LogInfo("AlcaBeamSpotManager")
00276 << "close payload because of movement in"
00277 << " X=" << deltaX
00278 << ", Y=" << deltaY
00279 << ", Z=" << deltaZ
00280 << ", sigmaZ=" << deltasigmaZ
00281 << ", dxdz=" << deltadxdz
00282 << ", dydz=" << deltadydz
00283 << ", widthX=" << deltawidthX
00284 << ", widthY=" << deltawidthY
00285 << endl;
00286 }
00287
00288 }
00289 if(docreate){
00290 if(foundShift){
00291 tmpBeamSpotMap_[firstBS->first] = weight(firstBS,nextBS);
00292 if (endOfRun){
00293
00294
00295 tmpBeamSpotMap_[nextBS->first] = nextBS->second;
00296 }
00297 }
00298 else if(!foundShift && !endOfRun){
00299 tmpBeamSpotMap_[firstBS->first] = weight(firstBS,nextBS);
00300 }
00301 else {
00302 tmpBeamSpotMap_[firstBS->first] = weight(firstBS,beamSpotMap_.end());
00303 }
00304 firstBS = nextBS;
00305 countlumi = 0;
00306
00307 }
00308
00309 ++countlumi;
00310
00311 currentBS = nextBS;
00312 nextBS = nextNextBS;
00313 nextNextBS++;
00314 ++iteration;
00315 }
00316 beamSpotMap_.clear();
00317 beamSpotMap_ = tmpBeamSpotMap_;
00318 }
00319 else if(beamSpotOutputBase_ == "runbased"){
00320 BeamSpot aBeamSpot = weight(beamSpotMap_.begin(),beamSpotMap_.end());
00321 LuminosityBlockNumber_t firstLumi = beamSpotMap_.begin()->first;
00322 beamSpotMap_.clear();
00323 beamSpotMap_[firstLumi] = aBeamSpot;
00324 }
00325 else{
00326 LogInfo("AlcaBeamSpotManager")
00327 << "Unrecognized BeamSpotOutputBase parameter: " << beamSpotOutputBase_
00328 << endl;
00329 }
00330 }
00331
00332
00333 BeamSpot AlcaBeamSpotManager::weight(const bsMap_iterator& begin,
00334 const bsMap_iterator& end){
00335 double x,xError = 0;
00336 double y,yError = 0;
00337 double z,zError = 0;
00338 double sigmaZ,sigmaZError = 0;
00339 double dxdz,dxdzError = 0;
00340 double dydz,dydzError = 0;
00341 double widthX,widthXError = 0;
00342 double widthY,widthYError = 0;
00343 LogInfo("AlcaBeamSpotManager")
00344 << "Weighted BeamSpot will span lumi "
00345 << begin->first << " to " << end->first
00346 << endl;
00347
00348 BeamSpot::BeamType type = BeamSpot::Unknown;
00349 for(bsMap_iterator it=begin; it!=end; it++){
00350 weight(x , xError , it->second.x0() , it->second.x0Error());
00351 weight(y , yError , it->second.y0() , it->second.y0Error());
00352 weight(z , zError , it->second.z0() , it->second.z0Error());
00353 weight(sigmaZ, sigmaZError, it->second.sigmaZ() , it->second.sigmaZ0Error());
00354 weight(dxdz , dxdzError , it->second.dxdz() , it->second.dxdzError());
00355 weight(dydz , dydzError , it->second.dydz() , it->second.dydzError());
00356 weight(widthX, widthXError, it->second.BeamWidthX(), it->second.BeamWidthXError());
00357 weight(widthY, widthYError, it->second.BeamWidthY(), it->second.BeamWidthYError());
00358 if(it->second.type() == BeamSpot::Tracker){
00359 type = BeamSpot::Tracker;
00360 }
00361 }
00362 BeamSpot::Point bsPosition(x,y,z);
00363 BeamSpot::CovarianceMatrix error;
00364 error(0,0) = xError*xError;
00365 error(1,1) = yError*yError;
00366 error(2,2) = zError*zError;
00367 error(3,3) = sigmaZError*sigmaZError;
00368 error(4,4) = dxdzError*dxdzError;
00369 error(5,5) = dydzError*dydzError;
00370 error(6,6) = widthXError*widthXError;
00371 BeamSpot weightedBeamSpot(bsPosition,sigmaZ,dxdz,dydz,widthX,error,type);
00372 weightedBeamSpot.setBeamWidthY(widthY);
00373 LogInfo("AlcaBeamSpotManager")
00374 << "Weighted BeamSpot will be:" <<'\n'
00375 << weightedBeamSpot
00376 << endl;
00377 return weightedBeamSpot;
00378 }
00379
00380
00381 void AlcaBeamSpotManager::weight(double& mean,double& meanError,const double& val,const double& valError){
00382 double tmpError = 0;
00383 if (meanError < 1e-8){
00384 tmpError = 1/(valError*valError);
00385 mean = val*tmpError;
00386 }
00387 else{
00388 tmpError = 1/(meanError*meanError) + 1/(valError*valError);
00389 mean = mean/(meanError*meanError) + val/(valError*valError);
00390 }
00391 mean = mean/tmpError;
00392 meanError = sqrt(1/tmpError);
00393 }
00394
00395
00396 pair<float,float> AlcaBeamSpotManager::delta(const float& x, const float& xError, const float& nextX, const float& nextXError){
00397 return pair<float,float>(x - nextX, sqrt(pow(xError,2) + pow(nextXError,2)) );
00398 }
00399
00400
00401 float AlcaBeamSpotManager::deltaSig(const float& num, const float& den){
00402 if(den != 0){
00403 return fabs(num/den);
00404 }
00405 else{
00406 return float(LONG_MAX);
00407 }
00408 }
00409
1.7.3