00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "DQM/BeamMonitor/plugins/BeamMonitorBx.h"
00011 #include "FWCore/ServiceRegistry/interface/Service.h"
00012 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00013 #include "DataFormats/Common/interface/View.h"
00014 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
00015 #include "FWCore/Framework/interface/MakerMacros.h"
00016 #include "FWCore/Framework/interface/Run.h"
00017 #include "FWCore/Framework/interface/LuminosityBlock.h"
00018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00019 #include <numeric>
00020 #include <math.h>
00021 #include <TMath.h>
00022 #include <iostream>
00023 #include <TStyle.h>
00024
00025 using namespace std;
00026 using namespace edm;
00027 using namespace reco;
00028
00029 const char * BeamMonitorBx::formatFitTime( const time_t & t ) {
00030 #define CET (+1)
00031 #define CEST (+2)
00032
00033 static char ts[] = "yyyy-Mm-dd hh:mm:ss";
00034 tm * ptm;
00035 ptm = gmtime ( &t );
00036 sprintf( ts, "%4d-%02d-%02d %02d:%02d:%02d", ptm->tm_year,ptm->tm_mon+1,ptm->tm_mday,(ptm->tm_hour+CEST)%24, ptm->tm_min, ptm->tm_sec);
00037
00038 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
00039 unsigned int b = strlen(ts);
00040 while (ts[--b] == ' ') {ts[b] = 0;}
00041 #endif
00042 return ts;
00043
00044 }
00045
00046
00047
00048
00049 BeamMonitorBx::BeamMonitorBx( const ParameterSet& ps ) :
00050 countBx_(0),countEvt_(0),countLumi_(0),resetHistos_(false) {
00051
00052 parameters_ = ps;
00053 monitorName_ = parameters_.getUntrackedParameter<string>("monitorName","YourSubsystemName");
00054 bsSrc_ = parameters_.getUntrackedParameter<InputTag>("beamSpot");
00055 fitNLumi_ = parameters_.getUntrackedParameter<int>("fitEveryNLumi",-1);
00056 resetFitNLumi_ = parameters_.getUntrackedParameter<int>("resetEveryNLumi",-1);
00057
00058 dbe_ = Service<DQMStore>().operator->();
00059
00060 if (monitorName_ != "" ) monitorName_ = monitorName_+"/" ;
00061
00062 theBeamFitter = new BeamFitter(parameters_);
00063 theBeamFitter->resetTrkVector();
00064 theBeamFitter->resetLSRange();
00065 theBeamFitter->resetRefTime();
00066 theBeamFitter->resetPVFitter();
00067
00068 if (fitNLumi_ <= 0) fitNLumi_ = 1;
00069 beginLumiOfBSFit_ = endLumiOfBSFit_ = 0;
00070 refBStime[0] = refBStime[1] = 0;
00071 lastlumi_ = 0;
00072 nextlumi_ = 0;
00073 firstlumi_ = 0;
00074 processed_ = false;
00075 countGoodFit_ = 0;
00076 }
00077
00078
00079 BeamMonitorBx::~BeamMonitorBx() {
00080 delete theBeamFitter;
00081 }
00082
00083
00084
00085 void BeamMonitorBx::beginJob() {
00086 varMap["x0_bx"] = "X_{0}";
00087 varMap["y0_bx"] = "Y_{0}";
00088 varMap["z0_bx"] = "Z_{0}";
00089 varMap["sigmaX_bx"] = "#sigma_{X}";
00090 varMap["sigmaY_bx"] = "#sigma_{Y}";
00091 varMap["sigmaZ_bx"] = "#sigma_{Z}";
00092
00093 varMap1["x"] = "X_{0} [cm]";
00094 varMap1["y"] = "Y_{0} [cm]";
00095 varMap1["z"] = "Z_{0} [cm]";
00096 varMap1["sigmaX"] = "#sigma_{X} [cm]";
00097 varMap1["sigmaY"] = "#sigma_{Y} [cm]";
00098 varMap1["sigmaZ"] = "#sigma_{Z} [cm]";
00099
00100
00101 dbe_->setCurrentFolder(monitorName_+"FitBx");
00102
00103 BookTables(1,varMap,"");
00104
00105
00106
00107 for (std::map<std::string,std::string>::const_iterator varName = varMap1.begin();
00108 varName != varMap1.end(); ++varName) {
00109 string subDir_ = "FitBx";
00110 subDir_ += "/";
00111 subDir_ += "All_";
00112 subDir_ += varName->first;
00113 dbe_->setCurrentFolder(monitorName_+subDir_);
00114 }
00115 dbe_->setCurrentFolder(monitorName_+"FitBx/All_nPVs");
00116 }
00117
00118
00119 void BeamMonitorBx::beginRun(const edm::Run& r, const EventSetup& context) {
00120
00121 ftimestamp = r.beginTime().value();
00122 tmpTime = ftimestamp >> 32;
00123 startTime = refTime = tmpTime;
00124
00125 }
00126
00127
00128 void BeamMonitorBx::beginLuminosityBlock(const LuminosityBlock& lumiSeg,
00129 const EventSetup& context) {
00130 int nthlumi = lumiSeg.luminosityBlock();
00131 const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
00132 const std::time_t ftmptime = fbegintimestamp >> 32;
00133
00134 if (countLumi_ == 0) {
00135 beginLumiOfBSFit_ = nthlumi;
00136 refBStime[0] = ftmptime;
00137 }
00138 if (beginLumiOfBSFit_ == 0) beginLumiOfBSFit_ = nextlumi_;
00139
00140 if (nthlumi < nextlumi_) return;
00141
00142 if (nthlumi > nextlumi_) {
00143 if (countLumi_ != 0 && processed_) {
00144 FitAndFill(lumiSeg,lastlumi_,nextlumi_,nthlumi);
00145 }
00146 nextlumi_ = nthlumi;
00147 edm::LogInfo("LS|BX|BeamMonitorBx") << "Next Lumi to Fit: " << nextlumi_ << endl;
00148 if (refBStime[0] == 0) refBStime[0] = ftmptime;
00149 }
00150 countLumi_++;
00151 if (processed_) processed_ = false;
00152 edm::LogInfo("LS|BX|BeamMonitorBx") << "Begin of Lumi: " << nthlumi << endl;
00153 }
00154
00155
00156 void BeamMonitorBx::analyze(const Event& iEvent,
00157 const EventSetup& iSetup ) {
00158 bool readEvent_ = true;
00159 const int nthlumi = iEvent.luminosityBlock();
00160 if (nthlumi < nextlumi_) {
00161 readEvent_ = false;
00162 edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from previous lumi section!" << endl;
00163 }
00164 if (nthlumi > nextlumi_) {
00165 readEvent_ = false;
00166 edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from next lumi section!!!" << endl;
00167 }
00168
00169 if (readEvent_) {
00170 countEvt_++;
00171 theBeamFitter->readEvent(iEvent);
00172 processed_ = true;
00173 }
00174
00175 }
00176
00177
00178
00179 void BeamMonitorBx::endLuminosityBlock(const LuminosityBlock& lumiSeg,
00180 const EventSetup& iSetup) {
00181 int nthlumi = lumiSeg.id().luminosityBlock();
00182 edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the last event before endLuminosityBlock: " << nthlumi << endl;
00183
00184 if (nthlumi < nextlumi_) return;
00185 const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
00186 const std::time_t fendtime = fendtimestamp >> 32;
00187 tmpTime = refBStime[1] = fendtime;
00188 }
00189
00190
00191
00192 void BeamMonitorBx::BookTables(int nBx, map<string,string> & vMap, string suffix_) {
00193
00194 dbe_->cd(monitorName_+"FitBx");
00195
00196 for (std::map<std::string,std::string>::const_iterator varName = vMap.begin();
00197 varName != vMap.end(); ++varName) {
00198 string tmpName = varName->first;
00199 if (!suffix_.empty()) {
00200 tmpName += "_";
00201 tmpName += suffix_;
00202 }
00203 if (dbe_->get(monitorName_+"FitBx/"+tmpName)) {
00204 edm::LogInfo("BX|BeamMonitorBx") << "Rebinning " << tmpName << endl;
00205 dbe_->removeElement(tmpName);
00206 }
00207
00208 hs[tmpName]=dbe_->book2D(tmpName,varName->second,3,0,3,nBx,0,nBx);
00209 hs[tmpName]->setBinLabel(1,"bx",1);
00210 hs[tmpName]->setBinLabel(2,varName->second,1);
00211 hs[tmpName]->setBinLabel(3,"#Delta "+varName->second,1);
00212 for (int i=0;i<nBx;i++) {
00213 hs[tmpName]->setBinLabel(i+1," ",2);
00214 }
00215 hs[tmpName]->getTH1()->SetOption("text");
00216 hs[tmpName]->Reset();
00217 }
00218 }
00219
00220
00221 void BeamMonitorBx::BookTrendHistos(bool plotPV,int nBx,map<string,string> & vMap,
00222 string subDir_, TString prefix_, TString suffix_) {
00223 int nType_ = 2;
00224 if (plotPV) nType_ = 4;
00225 std::ostringstream ss;
00226 std::ostringstream ss1;
00227 ss << setfill ('0') << setw (5) << nBx;
00228 ss1 << nBx;
00229
00230 for (int i = 0; i < nType_; i++) {
00231 for (std::map<std::string,std::string>::const_iterator varName = vMap.begin();
00232 varName != vMap.end(); ++varName) {
00233 string tmpDir_ = subDir_ + "/All_" + varName->first;
00234 dbe_->cd(monitorName_+tmpDir_);
00235 TString histTitle(varName->first);
00236 string tmpName;
00237 if (prefix_ != "") tmpName = prefix_ + "_" + varName->first;
00238 if (suffix_ != "") tmpName = tmpName + "_" + suffix_;
00239 tmpName = tmpName + "_" + ss.str();
00240
00241 TString histName(tmpName);
00242 string ytitle(varName->second);
00243 string xtitle("");
00244 string options("E1");
00245 bool createHisto = true;
00246 switch(i) {
00247 case 1:
00248 histName.Insert(histName.Index("_bx_",4),"_time");
00249 xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
00250 if (ytitle.find("sigma") == string::npos)
00251 histTitle += " coordinate of beam spot vs time (Fit)";
00252 else
00253 histTitle = histTitle.Insert(5," ") + " of beam spot vs time (Fit)";
00254 break;
00255 case 2:
00256 if (ytitle.find("sigma") == string::npos) {
00257 histName.Insert(0,"PV");
00258 histName.Insert(histName.Index("_bx_",4),"_lumi");
00259 histTitle.Insert(0,"Avg. ");
00260 histTitle += " position of primary vtx vs lumi";
00261 xtitle = "Lumisection [Bx# " + ss1.str() + "]";
00262 ytitle.insert(0,"PV");
00263 ytitle += " #pm #sigma_{PV";
00264 ytitle += varName->first;
00265 ytitle += "} (cm)";
00266 }
00267 else createHisto = false;
00268 break;
00269 case 3:
00270 if (ytitle.find("sigma") == string::npos) {
00271 histName.Insert(0,"PV");
00272 histName.Insert(histName.Index("_bx_",4),"_time");
00273 histTitle.Insert(0,"Avg. ");
00274 histTitle += " position of primary vtx vs time";
00275 xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
00276 ytitle.insert(0,"PV");
00277 ytitle += " #pm #sigma_{PV";
00278 ytitle += varName->first;
00279 ytitle += "} (cm)";
00280 }
00281 else createHisto = false;
00282 break;
00283 default:
00284 histName.Insert(histName.Index("_bx_",4),"_lumi");
00285 xtitle = "Lumisection [Bx# " + ss1.str() + "]";
00286 if (ytitle.find("sigma") == string::npos)
00287 histTitle += " coordinate of beam spot vs lumi (Fit)";
00288 else
00289 histTitle = histTitle.Insert(5," ") + " of beam spot vs lumi (Fit)";
00290 break;
00291 }
00292
00293 if (dbe_->get(monitorName_+tmpDir_+"/"+string(histName))) continue;
00294
00295 if (createHisto) {
00296 edm::LogInfo("BX|BeamMonitorBx") << "histName = " << histName << "; histTitle = " << histTitle << std::endl;
00297 hst[histName] = dbe_->book1D(histName,histTitle,40,0.5,40.5);
00298 hst[histName]->getTH1()->SetBit(TH1::kCanRebin);
00299 hst[histName]->setAxisTitle(xtitle,1);
00300 hst[histName]->setAxisTitle(ytitle,2);
00301 hst[histName]->getTH1()->SetOption("E1");
00302 if (histName.Contains("time")) {
00303 hst[histName]->getTH1()->SetBins(3600,0.5,3600+0.5);
00304 hst[histName]->setAxisTimeDisplay(1);
00305 hst[histName]->setAxisTimeFormat("%H:%M:%S",1);
00306 const char* eventTime = formatFitTime(startTime);
00307 TDatime da(eventTime);
00308 if (debug_) {
00309 edm::LogInfo("BX|BeamMonitorBx") << "TimeOffset = ";
00310 da.Print();
00311 }
00312 hst[histName]->getTH1()->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
00313 }
00314 }
00315 }
00316 }
00317
00318
00319 dbe_->cd(monitorName_+subDir_ + "/All_nPVs");
00320 TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
00321 string xtitle = "Lumisection [Bx# " + ss1.str() + "]";
00322
00323 hst[histName] = dbe_->book1D(histName,"Number of Good Reconstructed Vertices",40,0.5,40.5);
00324 hst[histName]->setAxisTitle(xtitle,1);
00325 hst[histName]->getTH1()->SetBit(TH1::kCanRebin);
00326 hst[histName]->getTH1()->SetOption("E1");
00327
00328 }
00329
00330
00331 void BeamMonitorBx::FitAndFill(const LuminosityBlock& lumiSeg,
00332 int &lastlumi,int &nextlumi,int &nthlumi){
00333 if (nthlumi <= nextlumi) return;
00334
00335 int currentlumi = nextlumi;
00336 edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the current fit: " << currentlumi << endl;
00337 lastlumi = currentlumi;
00338 endLumiOfBSFit_ = currentlumi;
00339
00340 edm::LogInfo("BX|BeamMonitorBx") << "Time lapsed = " << tmpTime - refTime << std:: endl;
00341
00342 if (resetHistos_) {
00343 edm::LogInfo("BX|BeamMonitorBx") << "Resetting Histograms" << endl;
00344 theBeamFitter->resetCutFlow();
00345 resetHistos_ = false;
00346 }
00347
00348 if (fitNLumi_ > 0)
00349 if (currentlumi%fitNLumi_!=0) return;
00350
00351 std::pair<int,int> fitLS = theBeamFitter->getFitLSRange();
00352 edm::LogInfo("LS|BX|BeamMonitorBx") << " [Fitter] Do BeamSpot Fit for LS = " << fitLS.first
00353 << " to " << fitLS.second << endl;
00354 edm::LogInfo("LS|BX|BeamMonitorBx") << " [BX] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_
00355 << " to " << endLumiOfBSFit_ << endl;
00356
00357 edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[0] = " << refBStime[0]
00358 << "; address = " << &refBStime[0] << std::endl;
00359 edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[1] = " << refBStime[1]
00360 << "; address = " << &refBStime[1] << std::endl;
00361
00362 if (theBeamFitter->runPVandTrkFitter()) {
00363 countGoodFit_++;
00364 edm::LogInfo("BX|BeamMonitorBx") << "Number of good fit = " << countGoodFit_ << endl;
00365 BeamSpotMapBx bsmap = theBeamFitter->getBeamSpotMap();
00366 std::map<int,int> npvsmap = theBeamFitter->getNPVsperBX();
00367 edm::LogInfo("BX|BeamMonitorBx") << "Number of bx = " << bsmap.size() << endl;
00368 if (bsmap.size() == 0) return;
00369 if (countBx_ < bsmap.size()) {
00370 countBx_ = bsmap.size();
00371 BookTables(countBx_,varMap,"");
00372 BookTables(countBx_,varMap,"all");
00373 for (BeamSpotMapBx::const_iterator abspot = bsmap.begin();
00374 abspot!= bsmap.end(); ++abspot) {
00375 int bx = abspot->first;
00376 BookTrendHistos(false,bx,varMap1,"FitBx","Trending","bx");
00377 }
00378 }
00379
00380 std::pair<int,int> LSRange = theBeamFitter->getFitLSRange();
00381 char tmpTitle[50];
00382 sprintf(tmpTitle,"%s %i %s %i %s"," [cm] (LS: ",LSRange.first," to ",LSRange.second,")");
00383 for (std::map<std::string,std::string>::const_iterator varName = varMap.begin();
00384 varName != varMap.end(); ++varName) {
00385 hs[varName->first]->setTitle(varName->second + " " + tmpTitle);
00386 hs[varName->first]->Reset();
00387 }
00388
00389 if (countGoodFit_ == 1)
00390 firstlumi_ = LSRange.first;
00391
00392 if (resetFitNLumi_ > 0 ) {
00393 char tmpTitle1[50];
00394 if ( countGoodFit_ > 1)
00395 sprintf(tmpTitle1,"%s %i %s %i %s"," [cm] (LS: ",firstlumi_," to ",LSRange.second,") [weighted average]");
00396 else
00397 sprintf(tmpTitle1,"%s","Need at least two fits to calculate weighted average");
00398 for (std::map<std::string,std::string>::const_iterator varName = varMap.begin();
00399 varName != varMap.end(); ++varName) {
00400 TString tmpName = varName->first + "_all";
00401 hs[tmpName]->setTitle(varName->second + " " + tmpTitle1);
00402 hs[tmpName]->Reset();
00403 }
00404 }
00405
00406 int nthBin = countBx_;
00407 for (BeamSpotMapBx::const_iterator abspot = bsmap.begin();
00408 abspot!= bsmap.end(); ++abspot,nthBin--) {
00409 reco::BeamSpot bs = abspot->second;
00410 int bx = abspot->first;
00411 int nPVs = npvsmap.find(bx)->second;
00412 edm::LogInfo("BeamMonitorBx") << "Number of PVs of bx#[" << bx << "] = " << nPVs << endl;
00413
00414 FillTables(bx,nthBin,varMap,bs,"");
00415
00416 FillTrendHistos(bx,nPVs,varMap1,bs,"Trending");
00417 }
00418
00419 weight(fbspotMap,bsmap);
00420
00421 nthBin = countBx_;
00422 if (resetFitNLumi_ > 0 && countGoodFit_ > 1) {
00423 for (BeamSpotMapBx::const_iterator abspot = fbspotMap.begin();
00424 abspot!= fbspotMap.end(); ++abspot,nthBin--) {
00425 reco::BeamSpot bs = abspot->second;
00426 int bx = abspot->first;
00427 FillTables(bx,nthBin,varMap,bs,"all");
00428 }
00429 }
00430 }
00431
00432
00433
00434 if (resetFitNLumi_ > 0 && currentlumi%resetFitNLumi_ == 0) {
00435 edm::LogInfo("LS|BX|BeamMonitorBx") << "Reset track collection for beam fit!!!" <<endl;
00436 resetHistos_ = true;
00437 theBeamFitter->resetTrkVector();
00438 theBeamFitter->resetLSRange();
00439 theBeamFitter->resetRefTime();
00440 theBeamFitter->resetPVFitter();
00441 beginLumiOfBSFit_ = 0;
00442 refBStime[0] = 0;
00443 }
00444 }
00445
00446
00447 void BeamMonitorBx::FillTables(int nthbx,int nthbin_,
00448 map<string,string>&vMap,
00449 reco::BeamSpot&bs_, string suffix_){
00450 map<string,pair<double,double> > val_;
00451 val_["x0_bx"] = pair<double,double>(bs_.x0(),bs_.x0Error());
00452 val_["y0_bx"] = pair<double,double>(bs_.y0(),bs_.y0Error());
00453 val_["z0_bx"] = pair<double,double>(bs_.z0(),bs_.z0Error());
00454 val_["sigmaX_bx"] = pair<double,double>(bs_.BeamWidthX(),bs_.BeamWidthXError());
00455 val_["sigmaY_bx"] = pair<double,double>(bs_.BeamWidthY(),bs_.BeamWidthYError());
00456 val_["sigmaZ_bx"] = pair<double,double>(bs_.sigmaZ(),bs_.sigmaZ0Error());
00457
00458 for (std::map<std::string,std::string>::const_iterator varName = vMap.begin();
00459 varName != vMap.end(); ++varName) {
00460 TString tmpName = varName->first;
00461 if (suffix_ != "") tmpName += TString("_"+suffix_);
00462 hs[tmpName]->setBinContent(1,nthbin_,nthbx);
00463 hs[tmpName]->setBinContent(2,nthbin_,val_[varName->first].first);
00464 hs[tmpName]->setBinContent(3,nthbin_,val_[varName->first].second);
00465 }
00466 }
00467
00468
00469 void BeamMonitorBx::FillTrendHistos(int nthBx, int nPVs_, map<string,string> & vMap,
00470 reco::BeamSpot & bs_, TString prefix_) {
00471 map<TString,pair<double,double> > val_;
00472 val_[TString(prefix_+"_x")] = pair<double,double>(bs_.x0(),bs_.x0Error());
00473 val_[TString(prefix_+"_y")] = pair<double,double>(bs_.y0(),bs_.y0Error());
00474 val_[TString(prefix_+"_z")] = pair<double,double>(bs_.z0(),bs_.z0Error());
00475 val_[TString(prefix_+"_sigmaX")] = pair<double,double>(bs_.BeamWidthX(),bs_.BeamWidthXError());
00476 val_[TString(prefix_+"_sigmaY")] = pair<double,double>(bs_.BeamWidthY(),bs_.BeamWidthYError());
00477 val_[TString(prefix_+"_sigmaZ")] = pair<double,double>(bs_.sigmaZ(),bs_.sigmaZ0Error());
00478
00479 std::ostringstream ss;
00480 ss << setfill ('0') << setw (5) << nthBx;
00481 int ntbin_ = tmpTime - startTime;
00482 for (map<TString,MonitorElement*>::iterator itHst = hst.begin();
00483 itHst != hst.end(); ++itHst) {
00484 if (!(itHst->first.Contains(ss.str()))) continue;
00485 if (itHst->first.Contains("nPVs")) continue;
00486 edm::LogInfo("BX|BeamMonitorBx") << "Filling histogram: " << itHst->first << endl;
00487 if (itHst->first.Contains("time")) {
00488 int idx = itHst->first.Index("_time",5);
00489 itHst->second->setBinContent(ntbin_,val_[itHst->first(0,idx)].first);
00490 itHst->second->setBinError(ntbin_,val_[itHst->first(0,idx)].second);
00491 }
00492 if (itHst->first.Contains("lumi")) {
00493 int idx = itHst->first.Index("_lumi",5);
00494 itHst->second->setBinContent(endLumiOfBSFit_,val_[itHst->first(0,idx)].first);
00495 itHst->second->setBinError(endLumiOfBSFit_,val_[itHst->first(0,idx)].second);
00496 }
00497 }
00498 TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
00499 if (hst[histName])
00500 hst[histName]->setBinContent(endLumiOfBSFit_,nPVs_);
00501 }
00502
00503
00504 void BeamMonitorBx::weight(BeamSpotMapBx& weightedMap_, const BeamSpotMapBx& newMap_){
00505 for (BeamSpotMapBx::const_iterator it = newMap_.begin();
00506 it != newMap_.end(); ++it) {
00507 if (weightedMap_.find(it->first) == weightedMap_.end() || (it->second.type() != 2)) {
00508 weightedMap_[it->first] = it->second;
00509 continue;
00510 }
00511
00512 BeamSpot obs = weightedMap_[it->first];
00513 double val_[8] = { obs.x0(),obs.y0(),obs.z0(),obs.sigmaZ(),
00514 obs.dxdz(),obs.dydz(),obs.BeamWidthX(),obs.BeamWidthY()};
00515 double valErr_[8] = { obs.x0Error(),obs.y0Error(),obs.z0Error(),
00516 obs.sigmaZ0Error(),obs.dxdzError(),obs.dydzError(),
00517 obs.BeamWidthXError(),obs.BeamWidthYError()};
00518
00519 reco::BeamSpot::BeamType type = reco::BeamSpot::Unknown;
00520 weight(val_[0], valErr_[0], it->second.x0() , it->second.x0Error());
00521 weight(val_[1], valErr_[1], it->second.y0() , it->second.y0Error());
00522 weight(val_[2], valErr_[2], it->second.z0() , it->second.z0Error());
00523 weight(val_[3], valErr_[3], it->second.sigmaZ() , it->second.sigmaZ0Error());
00524 weight(val_[4], valErr_[4], it->second.dxdz() , it->second.dxdzError());
00525 weight(val_[5], valErr_[5], it->second.dydz() , it->second.dydzError());
00526 weight(val_[6], valErr_[6], it->second.BeamWidthX(), it->second.BeamWidthXError());
00527 weight(val_[7], valErr_[7], it->second.BeamWidthY(), it->second.BeamWidthYError());
00528 if(it->second.type() == reco::BeamSpot::Tracker){
00529 type = reco::BeamSpot::Tracker;
00530 }
00531
00532 reco::BeamSpot::Point bsPosition(val_[0],val_[1],val_[2]);
00533 reco::BeamSpot::CovarianceMatrix error;
00534 for (int i=0;i<7;++i) {
00535 error(i,i) = valErr_[i]*valErr_[i];
00536 }
00537 reco::BeamSpot weightedBeamSpot(bsPosition,val_[3],val_[4],val_[5],val_[6],error,type);
00538 weightedBeamSpot.setBeamWidthY(val_[7]);
00539 LogInfo("BX|BeamMonitorBx")
00540 << weightedBeamSpot
00541 << endl;
00542 weightedMap_.erase(it->first);
00543 weightedMap_[it->first] = weightedBeamSpot;
00544 }
00545 }
00546
00547
00548 void BeamMonitorBx::weight(double& mean,double& meanError,const double& val,const double& valError){
00549 double tmpError = 0;
00550 if (meanError < 1e-8){
00551 tmpError = 1/(valError*valError);
00552 mean = val*tmpError;
00553 }
00554 else{
00555 tmpError = 1/(meanError*meanError) + 1/(valError*valError);
00556 mean = mean/(meanError*meanError) + val/(valError*valError);
00557 }
00558 mean = mean/tmpError;
00559 meanError = std::sqrt(1/tmpError);
00560 }
00561
00562
00563 void BeamMonitorBx::endRun(const Run& r, const EventSetup& context){
00564
00565 }
00566
00567
00568 void BeamMonitorBx::endJob(const LuminosityBlock& lumiSeg,
00569 const EventSetup& iSetup){
00570 }
00571
00572 DEFINE_FWK_MODULE(BeamMonitorBx);
1.7.3