00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "DTNoiseCalibration.h"
00011
00012
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/Run.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018
00019
00020 #include "Geometry/DTGeometry/interface/DTLayer.h"
00021 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00022 #include "Geometry/DTGeometry/interface/DTTopology.h"
00023 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00024
00025
00026 #include "DataFormats/DTDigi/interface/DTDigi.h"
00027 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
00028 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
00029 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
00030 #include "CondFormats/DTObjects/interface/DTReadOutMapping.h"
00031
00032
00033 #include "CondFormats/DTObjects/interface/DTTtrig.h"
00034 #include "CondFormats/DataRecord/interface/DTTtrigRcd.h"
00035
00036 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
00037
00038 #include "TH2F.h"
00039 #include "TFile.h"
00040
00041 using namespace edm;
00042 using namespace std;
00043
00044 DTNoiseCalibration::DTNoiseCalibration(const edm::ParameterSet& pset):
00045 digiLabel_( pset.getParameter<InputTag>("digiLabel") ),
00046 useTimeWindow_( pset.getParameter<bool>("useTimeWindow") ),
00047 triggerWidth_( pset.getParameter<double>("triggerWidth") ),
00048 timeWindowOffset_( pset.getParameter<int>("timeWindowOffset") ),
00049 maximumNoiseRate_( pset.getParameter<double>("maximumNoiseRate") ),
00050 useAbsoluteRate_( pset.getParameter<bool>("useAbsoluteRate") ),
00051 readDB_(true), defaultTtrig_(0),
00052 dbLabel_( pset.getUntrackedParameter<string>("dbLabel", "") ),
00053
00054 wireIdWithHisto_( std::vector<DTWireId>() ),
00055 lumiMax_(3000)
00056 {
00057
00058
00059
00060
00061
00062
00063
00064
00065 if( pset.exists("defaultTtrig") ){
00066 readDB_ = false;
00067 defaultTtrig_ = pset.getParameter<int>("defaultTtrig");
00068 }
00069
00070 if( pset.exists("cellsWithHisto") ){
00071 vector<string> cellsWithHisto = pset.getParameter<vector<string> >("cellsWithHisto");
00072 for(vector<string>::const_iterator cell = cellsWithHisto.begin(); cell != cellsWithHisto.end(); ++cell){
00073
00074 if( (*cell) != "" && (*cell) != "None"){
00075 stringstream linestr;
00076 int wheel,station,sector,sl,layer,wire;
00077 linestr << (*cell);
00078 linestr >> wheel >> station >> sector >> sl >> layer >> wire;
00079 wireIdWithHisto_.push_back(DTWireId(wheel,station,sector,sl,layer,wire));
00080 }
00081 }
00082 }
00083
00084
00085 string rootFileName = pset.getUntrackedParameter<string>("rootFileName","noise.root");
00086 rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
00087 rootFile_->cd();
00088 }
00089
00090 void DTNoiseCalibration::beginJob() {
00091 LogVerbatim("Calibration") << "[DTNoiseCalibration]: Begin job";
00092
00093 nevents_ = 0;
00094
00095 TH1::SetDefaultSumw2(true);
00096 int numBin = (triggerWidth_*32/25)/50;
00097 hTDCTriggerWidth_ = new TH1F("TDC_Time_Distribution", "TDC_Time_Distribution", numBin, 0, triggerWidth_*32/25);
00098 }
00099
00100 void DTNoiseCalibration::beginRun(const edm::Run& run, const edm::EventSetup& setup ) {
00101
00102
00103 setup.get<MuonGeometryRecord>().get(dtGeom_);
00104
00105
00106 if( readDB_ ) setup.get<DTTtrigRcd>().get(dbLabel_,tTrigMap_);
00107
00108 runBeginTime_ = time_t(run.beginTime().value()>>32);
00109 runEndTime_ = time_t(run.endTime().value()>>32);
00110
00111
00112
00113
00114
00115
00116
00117
00118 }
00119
00120 void DTNoiseCalibration::analyze(const edm::Event& event, const edm::EventSetup& setup){
00121 ++nevents_;
00122
00123
00124 Handle<DTDigiCollection> dtdigis;
00125 event.getByLabel(digiLabel_, dtdigis);
00126
00127
00128
00129
00130
00131
00132 time_t eventTime = time_t(event.time().value()>>32);
00133 unsigned int lumiSection = event.luminosityBlock();
00134
00135
00136 DTDigiCollection::DigiRangeIterator dtLayerId_It;
00137 for (dtLayerId_It=dtdigis->begin(); dtLayerId_It!=dtdigis->end(); ++dtLayerId_It){
00138
00139 float upperLimit = 0.;
00140 if(useTimeWindow_){
00141 if( readDB_ ){
00142 float tTrig,tTrigRMS,kFactor;
00143 DTSuperLayerId slId = ((*dtLayerId_It).first).superlayerId();
00144 int status = tTrigMap_->get( slId, tTrig, tTrigRMS, kFactor, DTTimeUnits::counts );
00145 if(status != 0) throw cms::Exception("DTNoiseCalibration") << "Could not find tTrig entry in DB for" << slId << endl;
00146 upperLimit = tTrig - timeWindowOffset_;
00147 } else {
00148 upperLimit = defaultTtrig_ - timeWindowOffset_;
00149 }
00150 }
00151
00152 double triggerWidth_s = 0.;
00153 if(useTimeWindow_) triggerWidth_s = ( (upperLimit*25)/32 )/1e9;
00154 else triggerWidth_s = double(triggerWidth_/1e9);
00155 LogTrace("Calibration") << ((*dtLayerId_It).first).superlayerId() << " Trigger width (s): " << triggerWidth_s;
00156
00157 for (DTDigiCollection::const_iterator digiIt = ((*dtLayerId_It).second).first;
00158 digiIt!=((*dtLayerId_It).second).second; ++digiIt){
00159
00160
00161 int tdcTime = (*digiIt).countsTDC();
00162 if( !useTimeWindow_ ){
00163 if( ( ((float)tdcTime*25)/32 ) > triggerWidth_ ){
00164 LogError("Calibration") << "Digi has a TDC time (ns) higher than the pre-defined TDC trigger width: " << ((float)tdcTime*25)/32;
00165 continue;
00166 }
00167 }
00168
00169 hTDCTriggerWidth_->Fill(tdcTime);
00170
00171 if( useTimeWindow_ && tdcTime > upperLimit) continue;
00172
00173
00174
00175
00176 const DTLayerId dtLId = (*dtLayerId_It).first;
00177 const DTTopology& dtTopo = dtGeom_->layer(dtLId)->specificTopology();
00178 const int firstWire = dtTopo.firstChannel();
00179 const int lastWire = dtTopo.lastChannel();
00180
00181 const int nWires = lastWire - firstWire + 1;
00182
00183
00184 if( theHistoOccupancyMap_.find(dtLId) == theHistoOccupancyMap_.end() ){
00185 string histoName = "DigiOccupancy_" + getLayerName(dtLId);
00186 rootFile_->cd();
00187 TH1F* hOccupancyHisto = new TH1F(histoName.c_str(), histoName.c_str(), nWires, firstWire, lastWire+1);
00188 LogTrace("Calibration") << " Created occupancy Histo: " << hOccupancyHisto->GetName();
00189 theHistoOccupancyMap_[dtLId] = hOccupancyHisto;
00190 }
00191 theHistoOccupancyMap_[dtLId]->Fill( (*digiIt).wire(), 1./triggerWidth_s );
00192
00193
00194 const DTChamberId dtChId = dtLId.chamberId();
00195 if( chamberOccupancyVsLumiMap_.find(dtChId) == chamberOccupancyVsLumiMap_.end() ){
00196 string histoName = "OccupancyVsLumi_" + getChamberName(dtChId);
00197 rootFile_->cd();
00198 TH1F* hOccupancyVsLumiHisto = new TH1F(histoName.c_str(), histoName.c_str(), lumiMax_, 0, lumiMax_);
00199 LogTrace("Calibration") << " Created occupancy histo: " << hOccupancyVsLumiHisto->GetName();
00200 chamberOccupancyVsLumiMap_[dtChId] = hOccupancyVsLumiHisto;
00201 }
00202 chamberOccupancyVsLumiMap_[dtChId]->Fill( lumiSection, 1./triggerWidth_s );
00203
00204 const DTWireId wireId(dtLId, (*digiIt).wire());
00205 if( find(wireIdWithHisto_.begin(),wireIdWithHisto_.end(),wireId) != wireIdWithHisto_.end() ){
00206 if( theHistoOccupancyVsLumiMap_.find(wireId) == theHistoOccupancyVsLumiMap_.end() ){
00207 string histoName = "OccupancyVsLumi_" + getChannelName(wireId);
00208 rootFile_->cd();
00209 TH1F* hOccupancyVsLumiHisto = new TH1F(histoName.c_str(), histoName.c_str(), lumiMax_, 0, lumiMax_);
00210 LogTrace("Calibration") << " Created occupancy histo: " << hOccupancyVsLumiHisto->GetName();
00211 theHistoOccupancyVsLumiMap_[wireId] = hOccupancyVsLumiHisto;
00212 }
00213 theHistoOccupancyVsLumiMap_[wireId]->Fill( lumiSection, 1./triggerWidth_s );
00214 }
00215
00216
00217 if( chamberOccupancyVsTimeMap_.find(dtChId) == chamberOccupancyVsTimeMap_.end() ){
00218 string histoName = "OccupancyVsTime_" + getChamberName(dtChId);
00219 float secPerBin = 20.0;
00220 unsigned int nBins = ( (unsigned int)(runEndTime_ - runBeginTime_) )/secPerBin;
00221 rootFile_->cd();
00222 TH1F* hOccupancyVsTimeHisto = new TH1F(histoName.c_str(), histoName.c_str(),
00223 nBins, (unsigned int)runBeginTime_,
00224 (unsigned int)runEndTime_);
00225 for(int k = 0; k < hOccupancyVsTimeHisto->GetNbinsX(); ++k){
00226 if( k%10 == 0 ){
00227 unsigned int binLowEdge = hOccupancyVsTimeHisto->GetBinLowEdge(k+1);
00228 time_t timeValue = time_t(binLowEdge);
00229 hOccupancyVsTimeHisto->GetXaxis()->SetBinLabel( (k+1),ctime(&timeValue) );
00230 }
00231 }
00232 size_t lastBin = hOccupancyVsTimeHisto->GetNbinsX();
00233 unsigned int binUpperEdge = hOccupancyVsTimeHisto->GetBinLowEdge(lastBin) +
00234 hOccupancyVsTimeHisto->GetBinWidth(lastBin);
00235 time_t timeValue = time_t(binUpperEdge);
00236 hOccupancyVsTimeHisto->GetXaxis()->SetBinLabel( (lastBin),ctime(&timeValue) );
00237
00238 LogTrace("Calibration") << " Created occupancy histo: " << hOccupancyVsTimeHisto->GetName();
00239 chamberOccupancyVsTimeMap_[dtChId] = hOccupancyVsTimeHisto;
00240 }
00241 chamberOccupancyVsTimeMap_[dtChId]->Fill( (unsigned int)eventTime, 1./triggerWidth_s );
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 }
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327 }
00328
00329 void DTNoiseCalibration::endJob(){
00330
00331
00332 LogVerbatim("Calibration") << "[DTNoiseCalibration] Total number of events analyzed: " << nevents_;
00333
00334
00335 rootFile_->cd();
00336 hTDCTriggerWidth_->Write();
00337
00338 double normalization = 1./double(nevents_);
00339
00340 for(map<DTWireId, TH1F*>::const_iterator wHisto = theHistoOccupancyVsLumiMap_.begin();
00341 wHisto != theHistoOccupancyVsLumiMap_.end(); ++wHisto){
00342 (*wHisto).second->Scale(normalization);
00343 (*wHisto).second->Write();
00344 }
00345
00346 for(map<DTChamberId, TH1F*>::const_iterator chHisto = chamberOccupancyVsLumiMap_.begin();
00347 chHisto != chamberOccupancyVsLumiMap_.end(); ++chHisto){
00348 (*chHisto).second->Scale(normalization);
00349 (*chHisto).second->Write();
00350 }
00351
00352 for(map<DTChamberId, TH1F*>::const_iterator chHisto = chamberOccupancyVsTimeMap_.begin();
00353 chHisto != chamberOccupancyVsTimeMap_.end(); ++chHisto){
00354 (*chHisto).second->Scale(normalization);
00355 (*chHisto).second->Write();
00356 }
00357
00358
00359 DTStatusFlag *statusMap = new DTStatusFlag();
00360 for(map<DTLayerId, TH1F*>::const_iterator lHisto = theHistoOccupancyMap_.begin();
00361 lHisto != theHistoOccupancyMap_.end();
00362 ++lHisto){
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 if((*lHisto).second){
00384 (*lHisto).second->Scale(normalization);
00385 rootFile_->cd();
00386 (*lHisto).second->Write();
00387 const DTTopology& dtTopo = dtGeom_->layer((*lHisto).first)->specificTopology();
00388 const int firstWire = dtTopo.firstChannel();
00389 const int lastWire = dtTopo.lastChannel();
00390
00391 const int nWires = lastWire - firstWire + 1;
00392
00393 double averageRate = 0.;
00394 for(int bin = 1; bin <= (*lHisto).second->GetNbinsX(); ++bin)
00395 averageRate += (*lHisto).second->GetBinContent(bin);
00396
00397 if(nWires) averageRate /= nWires;
00398 LogTrace("Calibration") << " Average rate = " << averageRate;
00399
00400 for(int i_wire = firstWire; i_wire <= lastWire; ++i_wire){
00401
00402 int bin = i_wire - firstWire + 1;
00403 double channelRate = (*lHisto).second->GetBinContent(bin);
00404 double rateOffset = (useAbsoluteRate_) ? 0. : averageRate;
00405 if( (channelRate - rateOffset) > maximumNoiseRate_ ){
00406 DTWireId wireID((*lHisto).first, i_wire);
00407 statusMap->setCellNoise(wireID,1);
00408 LogVerbatim("Calibration") << ">>> Channel noisy: " << wireID;
00409 }
00410 }
00411 }
00412 }
00413 LogVerbatim("Calibration") << "Writing noise map object to DB";
00414 string record = "DTStatusFlagRcd";
00415 DTCalibDBUtils::writeToDB<DTStatusFlag>(record, statusMap);
00416 }
00417
00418 DTNoiseCalibration::~DTNoiseCalibration(){
00419 rootFile_->Close();
00420 }
00421
00422 string DTNoiseCalibration::getChannelName(const DTWireId& wId) const{
00423 stringstream channelName;
00424 channelName << "Wh" << wId.wheel() << "_St" << wId.station() << "_Sec" << wId.sector()
00425 << "_SL" << wId.superlayer() << "_L" << wId.layer() << "_W"<< wId.wire();
00426
00427 return channelName.str();
00428 }
00429
00430 string DTNoiseCalibration::getLayerName(const DTLayerId& lId) const{
00431
00432 const DTSuperLayerId dtSLId = lId.superlayerId();
00433 const DTChamberId dtChId = dtSLId.chamberId();
00434 stringstream Layer; Layer << lId.layer();
00435 stringstream superLayer; superLayer << dtSLId.superlayer();
00436 stringstream wheel; wheel << dtChId.wheel();
00437 stringstream station; station << dtChId.station();
00438 stringstream sector; sector << dtChId.sector();
00439
00440 string layerName =
00441 "W" + wheel.str()
00442 + "_St" + station.str()
00443 + "_Sec" + sector.str()
00444 + "_SL" + superLayer.str()
00445 + "_L" + Layer.str();
00446
00447 return layerName;
00448 }
00449
00450 string DTNoiseCalibration::getSuperLayerName(const DTSuperLayerId& dtSLId) const{
00451
00452 const DTChamberId dtChId = dtSLId.chamberId();
00453 stringstream superLayer; superLayer << dtSLId.superlayer();
00454 stringstream wheel; wheel << dtChId.wheel();
00455 stringstream station; station << dtChId.station();
00456 stringstream sector; sector << dtChId.sector();
00457
00458 string superLayerName =
00459 "W" + wheel.str()
00460 + "_St" + station.str()
00461 + "_Sec" + sector.str()
00462 + "_SL" + superLayer.str();
00463
00464 return superLayerName;
00465 }
00466
00467 string DTNoiseCalibration::getChamberName(const DTChamberId& dtChId) const{
00468
00469 stringstream wheel; wheel << dtChId.wheel();
00470 stringstream station; station << dtChId.station();
00471 stringstream sector; sector << dtChId.sector();
00472
00473 string chamberName =
00474 "W" + wheel.str()
00475 + "_St" + station.str()
00476 + "_Sec" + sector.str();
00477
00478 return chamberName;
00479 }