DTCalibration/plugins/DTVDriftSegment.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author A. Vilela Pereira
0006  */
0007
0008 #include "DTVDriftSegment.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0016 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0017 #include "CondFormats/DTObjects/interface/DTMtime.h"
0018 #include "CondFormats/DataRecord/interface/DTMtimeRcd.h"
0019 #include "CondFormats/DTObjects/interface/DTRecoConditions.h"
0020 #include "CondFormats/DataRecord/interface/DTRecoConditionsVdriftRcd.h"
0021
0022 #include "CalibMuon/DTCalibration/interface/DTResidualFitter.h"
0023 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0024 #include "FWCore/Framework/interface/ConsumesCollector.h"
0025 #include <string>
0026 #include <vector>
0027
0028 #include "TH1F.h"
0029 #include "TFile.h"
0030
0031 using namespace std;
0032 using namespace edm;
0033
0034 namespace dtCalibration {
0035
0036   DTVDriftSegment::DTVDriftSegment(const ParameterSet& pset, edm::ConsumesCollector cc)
0037       : nSigmas_(pset.getUntrackedParameter<unsigned int>("nSigmasFitRange", 1)),
0038         mTimeMap_(nullptr),
0039         vDriftMap_(nullptr) {
0040     string rootFileName = pset.getParameter<string>("rootFileName");
0041     rootFile_ = new TFile(rootFileName.c_str(), "READ");
0042
0043     bool debug = pset.getUntrackedParameter<bool>("debug", false);
0044     fitter_ = new DTResidualFitter(debug);
0045     //if(debug) fitter_->setVerbosity(1);
0046
0047     readLegacyVDriftDB = pset.getParameter<bool>("readLegacyVDriftDB");
0048     if (readLegacyVDriftDB) {
0049       mTimeMapToken_ = cc.esConsumes<edm::Transition::BeginRun>();
0050     } else {
0051       vDriftMapToken_ = cc.esConsumes<edm::Transition::BeginRun>();
0052     }
0053   }
0054
0055   DTVDriftSegment::~DTVDriftSegment() {
0056     rootFile_->Close();
0057     delete fitter_;
0058   }
0059
0060   void DTVDriftSegment::setES(const edm::EventSetup& setup) {
0061     // Get the map of vdrift from the setup
0062     if (readLegacyVDriftDB) {
0063       mTimeMap_ = &setup.getData(mTimeMapToken_);
0064     } else {
0065       vDriftMap_ = &setup.getData(vDriftMapToken_);
0066       // Consistency check: no parametrization is implemented for the time being
0067       int version = vDriftMap_->version();
0068       if (version != 1) {
0069         throw cms::Exception("Configuration") << "only version 1 is presently supported for VDriftDB";
0070       }
0071     }
0072   }
0073
0074   DTVDriftData DTVDriftSegment::compute(DTSuperLayerId const& slId) {
0075     // Get original value from DB; vdrift is cm/ns , resolution is cm
0076     // Note that resolution is irrelevant as it is no longer used anywhere in reconstruction.
0077
0078     float vDrift = 0., resolution = 0.;
0079     if (readLegacyVDriftDB) {  // Legacy format
0080       int status = mTimeMap_->get(slId, vDrift, resolution, DTVelocityUnits::cm_per_ns);
0081       if (status != 0)
0082         throw cms::Exception("DTCalibration") << "Could not find vDrift entry in DB for" << slId << endl;
0083     } else {  // New DB format
0084       vDrift = vDriftMap_->get(DTWireId(slId.rawId()));
0085     }
0086
0087     // For RZ superlayers use original value
0088     if (slId.superLayer() == 2) {
0089       LogTrace("Calibration") << "[DTVDriftSegment]: RZ superlayer\n"
0090                               << "                   Will use original vDrift and resolution.";
0091       return DTVDriftData(vDrift, resolution);
0092     } else {
0093       TH1F* vDriftCorrHisto = getHisto(slId);
0094       // If empty histogram
0095       if (vDriftCorrHisto->GetEntries() == 0) {
0096         LogError("Calibration") << "[DTVDriftSegment]: Histogram " << vDriftCorrHisto->GetName() << " is empty.\n"
0097                                 << "                   Will use original vDrift and resolution.";
0098         return DTVDriftData(vDrift, resolution);
0099       }
0100
0101       LogTrace("Calibration") << "[DTVDriftSegment]: Fitting histogram " << vDriftCorrHisto->GetName();
0102       DTResidualFitResult fitResult = fitter_->fitResiduals(*vDriftCorrHisto, nSigmas_);
0103       LogTrace("Calibration") << "[DTVDriftSegment]: \n"
0104                               << "   Fit Mean  = " << fitResult.fitMean << " +/- " << fitResult.fitMeanError << "\n"
0105                               << "   Fit Sigma = " << fitResult.fitSigma << " +/- " << fitResult.fitSigmaError;
0106
0107       float vDriftCorr = fitResult.fitMean;
0108       float vDriftNew = vDrift * (1. - vDriftCorr);
0109       float resolutionNew = resolution;
0110       return DTVDriftData(vDriftNew, resolutionNew);
0111     }
0112   }
0113
0114   TH1F* DTVDriftSegment::getHisto(const DTSuperLayerId& slId) {
0115     string histoName = getHistoName(slId);
0116     TH1F* histo = static_cast<TH1F*>(rootFile_->Get(histoName.c_str()));
0117     if (!histo)
0118       throw cms::Exception("DTCalibration") << "v-drift correction histogram not found:" << histoName << endl;
0119     return histo;
0120   }
0121
0122   string DTVDriftSegment::getHistoName(const DTSuperLayerId& slId) {
0123     DTChamberId chId = slId.chamberId();
0124
0125     // Compose the chamber name
0126     std::string wheel = std::to_string(chId.wheel());
0127     std::string station = std::to_string(chId.station());
0128     std::string sector = std::to_string(chId.sector());
0129
0130     string chHistoName = "_W" + wheel + "_St" + station + "_Sec" + sector;
0131
0132     return (slId.superLayer() != 2) ? ("hRPhiVDriftCorr" + chHistoName) : ("hRZVDriftCorr" + chHistoName);
0133   }
0134
0135 }  // namespace dtCalibration