DTCalibration/plugins/DTNoiseCalibration.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Mila - INFN Torino
0005  *          A. Vilela Pereira - INFN Torino
0006  */
0007
0008 #include "DTNoiseCalibration.h"
0009
0010 // Framework
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Run.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 // Geometry
0018 #include "Geometry/DTGeometry/interface/DTLayer.h"
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 #include "Geometry/DTGeometry/interface/DTTopology.h"
0021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0022
0023 // Digis
0024 #include "DataFormats/DTDigi/interface/DTDigi.h"
0025 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0026 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0027 #include "CondFormats/DTObjects/interface/DTReadOutMapping.h"
0028
0029 // Database
0030 #include "CondFormats/DTObjects/interface/DTTtrig.h"
0031 #include "CondFormats/DataRecord/interface/DTTtrigRcd.h"
0032
0033 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0034
0035 #include "TH2F.h"
0036 #include "TFile.h"
0037
0038 using namespace edm;
0039 using namespace std;
0040
0041 DTNoiseCalibration::DTNoiseCalibration(const edm::ParameterSet& pset)
0042     : digiToken_(consumes<DTDigiCollection>(pset.getParameter<InputTag>("digiLabel"))),
0043       useTimeWindow_(pset.getParameter<bool>("useTimeWindow")),
0044       triggerWidth_(pset.getParameter<double>("triggerWidth")),
0045       timeWindowOffset_(pset.getParameter<int>("timeWindowOffset")),
0046       maximumNoiseRate_(pset.getParameter<double>("maximumNoiseRate")),
0047       useAbsoluteRate_(pset.getParameter<bool>("useAbsoluteRate")),
0048       readDB_(true),
0049       defaultTtrig_(0),
0050       //fastAnalysis_( pset.getParameter<bool>("fastAnalysis", true) ),
0051       wireIdWithHisto_(std::vector<DTWireId>()),
0052       lumiMax_(3000),
0053       dtGeomToken_(esConsumes<edm::Transition::BeginRun>()),
0054       ttrigToken_(
0055           esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", pset.getUntrackedParameter<string>("dbLabel")))) {
0056   // Get the debug parameter for verbose output
0057   //debug = ps.getUntrackedParameter<bool>("debug");
0058   /*// The analysis type
0059   // The wheel & sector interested for the time-dependent analysis
0060   wh = ps.getUntrackedParameter<int>("wheel", 0);
0061   sect = ps.getUntrackedParameter<int>("sector", 6);*/
0062
0063   if (pset.exists("defaultTtrig")) {
0064     readDB_ = false;
0065     defaultTtrig_ = pset.getParameter<int>("defaultTtrig");
0066   }
0067
0068   if (pset.exists("cellsWithHisto")) {
0069     vector<string> cellsWithHisto = pset.getParameter<vector<string> >("cellsWithHisto");
0070     for (vector<string>::const_iterator cell = cellsWithHisto.begin(); cell != cellsWithHisto.end(); ++cell) {
0071       //FIXME: Use regex to check whether format is right
0072       if ((!cell->empty()) && (*cell) != "None") {
0073         stringstream linestr;
0074         int wheel, station, sector, sl, layer, wire;
0075         linestr << (*cell);
0076         linestr >> wheel >> station >> sector >> sl >> layer >> wire;
0077         wireIdWithHisto_.push_back(DTWireId(wheel, station, sector, sl, layer, wire));
0078       }
0079     }
0080   }
0081
0082   // The root file which will contain the histos
0083   string rootFileName = pset.getUntrackedParameter<string>("rootFileName", "noise.root");
0084   rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
0085   rootFile_->cd();
0086 }
0087
0088 void DTNoiseCalibration::beginJob() {
0089   LogVerbatim("Calibration") << "[DTNoiseCalibration]: Begin job";
0090
0091   nevents_ = 0;
0092
0093   TH1::SetDefaultSumw2(true);
0094   int numBin = (triggerWidth_ * 32 / 25) / 50;
0095   hTDCTriggerWidth_ = new TH1F("TDC_Time_Distribution", "TDC_Time_Distribution", numBin, 0, triggerWidth_ * 32 / 25);
0096 }
0097
0098 void DTNoiseCalibration::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
0099   // Get the DT Geometry
0100   dtGeom_ = setup.getHandle(dtGeomToken_);
0101   // tTrig
0102   if (readDB_)
0103     tTrigMap_ = &setup.getData(ttrigToken_);
0104   runBeginTime_ = time_t(run.beginTime().value() >> 32);
0105   runEndTime_ = time_t(run.endTime().value() >> 32);
0106 }
0107
0108 void DTNoiseCalibration::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0109   ++nevents_;
0110
0111   // Get the digis from the event
0112   const Handle<DTDigiCollection>& dtdigis = event.getHandle(digiToken_);
0113
0114   time_t eventTime = time_t(event.time().value() >> 32);
0115   unsigned int lumiSection = event.luminosityBlock();
0116
0117   // Loop over digis
0118   DTDigiCollection::DigiRangeIterator dtLayerId_It;
0119   for (dtLayerId_It = dtdigis->begin(); dtLayerId_It != dtdigis->end(); ++dtLayerId_It) {
0120     // Define time window
0121     float upperLimit = 0.;
0122     if (useTimeWindow_) {
0123       if (readDB_) {
0124         float tTrig, tTrigRMS, kFactor;
0125         DTSuperLayerId slId = ((*dtLayerId_It).first).superlayerId();
0126         int status = tTrigMap_->get(slId, tTrig, tTrigRMS, kFactor, DTTimeUnits::counts);
0127         if (status != 0)
0128           throw cms::Exception("DTNoiseCalibration") << "Could not find tTrig entry in DB for" << slId << endl;
0129         upperLimit = tTrig - timeWindowOffset_;
0130       } else {
0131         upperLimit = defaultTtrig_ - timeWindowOffset_;
0132       }
0133     }
0134
0135     double triggerWidth_s = 0.;
0136     if (useTimeWindow_)
0137       triggerWidth_s = ((upperLimit * 25) / 32) / 1e9;
0138     else
0139       triggerWidth_s = double(triggerWidth_ / 1e9);
0140     LogTrace("Calibration") << ((*dtLayerId_It).first).superlayerId() << " Trigger width (s): " << triggerWidth_s;
0141
0142     for (DTDigiCollection::const_iterator digiIt = ((*dtLayerId_It).second).first;
0143          digiIt != ((*dtLayerId_It).second).second;
0144          ++digiIt) {
0145       //Check the TDC trigger width
0146       int tdcTime = (*digiIt).countsTDC();
0147       if (!useTimeWindow_) {
0148         if ((((float)tdcTime * 25) / 32) > triggerWidth_) {
0149           LogError("Calibration") << "Digi has a TDC time (ns) higher than the pre-defined TDC trigger width: "
0150                                   << ((float)tdcTime * 25) / 32;
0151           continue;
0152         }
0153       }
0154
0155       hTDCTriggerWidth_->Fill(tdcTime);
0156
0157       if (useTimeWindow_ && tdcTime > upperLimit)
0158         continue;
0159
0160       const DTLayerId dtLId = (*dtLayerId_It).first;
0161       const DTTopology& dtTopo = dtGeom_->layer(dtLId)->specificTopology();
0162       const int firstWire = dtTopo.firstChannel();
0163       const int lastWire = dtTopo.lastChannel();
0164       const int nWires = lastWire - firstWire + 1;
0165
0166       // Book the occupancy histos
0167       if (theHistoOccupancyMap_.find(dtLId) == theHistoOccupancyMap_.end()) {
0168         string histoName = "DigiOccupancy_" + getLayerName(dtLId);
0169         rootFile_->cd();
0170         TH1F* hOccupancyHisto = new TH1F(histoName.c_str(), histoName.c_str(), nWires, firstWire, lastWire + 1);
0171         LogTrace("Calibration") << "  Created occupancy Histo: " << hOccupancyHisto->GetName();
0172         theHistoOccupancyMap_[dtLId] = hOccupancyHisto;
0173       }
0174       theHistoOccupancyMap_[dtLId]->Fill((*digiIt).wire(), 1. / triggerWidth_s);
0175
0176       // Histos vs lumi section
0177       const DTChamberId dtChId = dtLId.chamberId();
0178       if (chamberOccupancyVsLumiMap_.find(dtChId) == chamberOccupancyVsLumiMap_.end()) {
0179         string histoName = "OccupancyVsLumi_" + getChamberName(dtChId);
0180         rootFile_->cd();
0181         TH1F* hOccupancyVsLumiHisto = new TH1F(histoName.c_str(), histoName.c_str(), lumiMax_, 0, lumiMax_);
0182         LogTrace("Calibration") << "  Created occupancy histo: " << hOccupancyVsLumiHisto->GetName();
0183         chamberOccupancyVsLumiMap_[dtChId] = hOccupancyVsLumiHisto;
0184       }
0185       chamberOccupancyVsLumiMap_[dtChId]->Fill(lumiSection, 1. / triggerWidth_s);
0186
0187       const DTWireId wireId(dtLId, (*digiIt).wire());
0188       if (find(wireIdWithHisto_.begin(), wireIdWithHisto_.end(), wireId) != wireIdWithHisto_.end()) {
0189         if (theHistoOccupancyVsLumiMap_.find(wireId) == theHistoOccupancyVsLumiMap_.end()) {
0190           string histoName = "OccupancyVsLumi_" + getChannelName(wireId);
0191           rootFile_->cd();
0192           TH1F* hOccupancyVsLumiHisto = new TH1F(histoName.c_str(), histoName.c_str(), lumiMax_, 0, lumiMax_);
0193           LogTrace("Calibration") << "  Created occupancy histo: " << hOccupancyVsLumiHisto->GetName();
0194           theHistoOccupancyVsLumiMap_[wireId] = hOccupancyVsLumiHisto;
0195         }
0196         theHistoOccupancyVsLumiMap_[wireId]->Fill(lumiSection, 1. / triggerWidth_s);
0197       }
0198
0199       // Histos vs time
0200       if (chamberOccupancyVsTimeMap_.find(dtChId) == chamberOccupancyVsTimeMap_.end()) {
0201         string histoName = "OccupancyVsTime_" + getChamberName(dtChId);
0202         float secPerBin = 20.0;
0203         unsigned int nBins = ((unsigned int)(runEndTime_ - runBeginTime_)) / secPerBin;
0204         rootFile_->cd();
0205         TH1F* hOccupancyVsTimeHisto = new TH1F(
0206             histoName.c_str(), histoName.c_str(), nBins, (unsigned int)runBeginTime_, (unsigned int)runEndTime_);
0207         for (int k = 0; k < hOccupancyVsTimeHisto->GetNbinsX(); ++k) {
0208           if (k % 10 == 0) {
0209             unsigned int binLowEdge = hOccupancyVsTimeHisto->GetBinLowEdge(k + 1);
0210             time_t timeValue = time_t(binLowEdge);
0211             hOccupancyVsTimeHisto->GetXaxis()->SetBinLabel((k + 1), ctime(&timeValue));
0212           }
0213         }
0214         size_t lastBin = hOccupancyVsTimeHisto->GetNbinsX();
0215         unsigned int binUpperEdge =
0216             hOccupancyVsTimeHisto->GetBinLowEdge(lastBin) + hOccupancyVsTimeHisto->GetBinWidth(lastBin);
0217         time_t timeValue = time_t(binUpperEdge);
0218         hOccupancyVsTimeHisto->GetXaxis()->SetBinLabel((lastBin), ctime(&timeValue));
0219
0220         LogTrace("Calibration") << "  Created occupancy histo: " << hOccupancyVsTimeHisto->GetName();
0221         chamberOccupancyVsTimeMap_[dtChId] = hOccupancyVsTimeHisto;
0222       }
0223       chamberOccupancyVsTimeMap_[dtChId]->Fill((unsigned int)eventTime, 1. / triggerWidth_s);
0224     }
0225   }
0226 }
0227
0228 void DTNoiseCalibration::endJob() {
0229   LogVerbatim("Calibration") << "[DTNoiseCalibration] Total number of events analyzed: " << nevents_;
0230
0231   // Save the TDC digi plot
0232   rootFile_->cd();
0233   hTDCTriggerWidth_->Write();
0234
0235   double normalization = 1. / double(nevents_);
0236
0237   for (map<DTWireId, TH1F*>::const_iterator wHisto = theHistoOccupancyVsLumiMap_.begin();
0238        wHisto != theHistoOccupancyVsLumiMap_.end();
0239        ++wHisto) {
0240     (*wHisto).second->Scale(normalization);
0241     (*wHisto).second->Write();
0242   }
0243
0244   for (map<DTChamberId, TH1F*>::const_iterator chHisto = chamberOccupancyVsLumiMap_.begin();
0245        chHisto != chamberOccupancyVsLumiMap_.end();
0246        ++chHisto) {
0247     (*chHisto).second->Scale(normalization);
0248     (*chHisto).second->Write();
0249   }
0250
0251   for (map<DTChamberId, TH1F*>::const_iterator chHisto = chamberOccupancyVsTimeMap_.begin();
0252        chHisto != chamberOccupancyVsTimeMap_.end();
0253        ++chHisto) {
0254     (*chHisto).second->Scale(normalization);
0255     (*chHisto).second->Write();
0256   }
0257
0258   // Save on file the occupancy histos and write the list of noisy cells
0259   DTStatusFlag statusMap;
0260   for (map<DTLayerId, TH1F*>::const_iterator lHisto = theHistoOccupancyMap_.begin();
0261        lHisto != theHistoOccupancyMap_.end();
0262        ++lHisto) {
0263     if ((*lHisto).second) {
0264       (*lHisto).second->Scale(normalization);
0265       rootFile_->cd();
0266       (*lHisto).second->Write();
0267       const DTTopology& dtTopo = dtGeom_->layer((*lHisto).first)->specificTopology();
0268       const int firstWire = dtTopo.firstChannel();
0269       const int lastWire = dtTopo.lastChannel();
0270       const int nWires = lastWire - firstWire + 1;
0271       // Find average in layer
0272       double averageRate = 0.;
0273       for (int bin = 1; bin <= (*lHisto).second->GetNbinsX(); ++bin)
0274         averageRate += (*lHisto).second->GetBinContent(bin);
0275
0276       if (nWires)
0277         averageRate /= nWires;
0278       LogTrace("Calibration") << "  Average rate = " << averageRate;
0279
0280       for (int i_wire = firstWire; i_wire <= lastWire; ++i_wire) {
0281         // From definition of "noisy cell"
0282         int bin = i_wire - firstWire + 1;
0283         double channelRate = (*lHisto).second->GetBinContent(bin);
0284         double rateOffset = (useAbsoluteRate_) ? 0. : averageRate;
0285         if ((channelRate - rateOffset) > maximumNoiseRate_) {
0286           DTWireId wireID((*lHisto).first, i_wire);
0287           statusMap.setCellNoise(wireID, true);
0288           LogVerbatim("Calibration") << ">>> Channel noisy: " << wireID;
0289         }
0290       }
0291     }
0292   }
0293   LogVerbatim("Calibration") << "Writing noise map object to DB";
0294   string record = "DTStatusFlagRcd";
0295   DTCalibDBUtils::writeToDB<DTStatusFlag>(record, statusMap);
0296 }
0297
0298 DTNoiseCalibration::~DTNoiseCalibration() { rootFile_->Close(); }
0299
0300 string DTNoiseCalibration::getChannelName(const DTWireId& wId) const {
0301   string channelName = "Wh" + std::to_string(wId.wheel()) + "_St" + std::to_string(wId.station()) + "_Sec" +
0302                        std::to_string(wId.sector()) + "_SL" + std::to_string(wId.superlayer()) + "_L" +
0303                        std::to_string(wId.layer()) + "_W" + std::to_string(wId.wire());
0304   return channelName;
0305 }
0306
0307 string DTNoiseCalibration::getLayerName(const DTLayerId& lId) const {
0308   const DTSuperLayerId dtSLId = lId.superlayerId();
0309   const DTChamberId dtChId = dtSLId.chamberId();
0310   string layerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0311                      std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer()) + "_L" +
0312                      std::to_string(lId.layer());
0313
0314   return layerName;
0315 }
0316
0317 string DTNoiseCalibration::getSuperLayerName(const DTSuperLayerId& dtSLId) const {
0318   const DTChamberId dtChId = dtSLId.chamberId();
0319
0320   string superLayerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0321                           std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer());
0322
0323   return superLayerName;
0324 }
0325
0326 string DTNoiseCalibration::getChamberName(const DTChamberId& dtChId) const {
0327   string chamberName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0328                        std::to_string(dtChId.sector());
0329
0330   return chamberName;
0331 }