DTCalibration/plugins/DTTTrigResidualCorrection.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author A. Vilela Pereira
0005  */
0006
0007 #include "DTTTrigResidualCorrection.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "DataFormats/MuonDetId/interface/DTSuperLayerId.h"
0014 #include "CondFormats/DTObjects/interface/DTTtrig.h"
0015 #include "CondFormats/DataRecord/interface/DTTtrigRcd.h"
0016 #include "CondFormats/DTObjects/interface/DTMtime.h"
0017 #include "CondFormats/DataRecord/interface/DTMtimeRcd.h"
0018 #include "CondFormats/DTObjects/interface/DTRecoConditions.h"
0019 #include "CondFormats/DataRecord/interface/DTRecoConditionsVdriftRcd.h"
0020
0021 #include "CalibMuon/DTCalibration/interface/DTResidualFitter.h"
0022
0023 #include "TFile.h"
0024 #include "TH1F.h"
0025 #include "TF1.h"
0026 #include "TCanvas.h"
0027
0028 #include "RooPlot.h"
0029 #include "RooRealVar.h"
0030 #include "RooDataHist.h"
0031 #include "RooGaussian.h"
0032 #include "RooAddPdf.h"
0033 #include "RooFitResult.h"
0034 #include "RooGlobalFunc.h"
0035
0036 #include <string>
0037 #include <sstream>
0038 #include <fstream>
0039
0040 using namespace std;
0041 using namespace edm;
0042
0043 namespace dtCalibration {
0044
0045   DTTTrigResidualCorrection::DTTTrigResidualCorrection(const ParameterSet& pset, edm::ConsumesCollector cc) {
0046     string residualsRootFile = pset.getParameter<string>("residualsRootFile");
0047     rootFile_ = new TFile(residualsRootFile.c_str(), "READ");
0048     rootBaseDir_ = pset.getUntrackedParameter<string>("rootBaseDir", "/DQMData/DT/DTCalibValidation");
0049     useFit_ = pset.getParameter<bool>("useFitToResiduals");
0050     //useConstantvDrift_ = pset.getParameter<bool>("useConstantDriftVelocity");
0051     useSlopesCalib_ = pset.getUntrackedParameter<bool>("useSlopesCalib", false);
0052     readLegacyVDriftDB = pset.getParameter<bool>("readLegacyVDriftDB");
0053     ttrigToken_ =
0054         cc.esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", pset.getUntrackedParameter<string>("dbLabel")));
0055     mTimeMapToken_ = cc.esConsumes<edm::Transition::BeginRun>();
0056     vDriftMapToken_ = cc.esConsumes<edm::Transition::BeginRun>();
0057
0058     // Load external slopes
0059     if (useSlopesCalib_) {
0060       ifstream fileSlopes;
0061       fileSlopes.open(pset.getParameter<FileInPath>("slopesFileName").fullPath().c_str());
0062
0063       int tmp_wheel = 0, tmp_sector = 0, tmp_station = 0, tmp_SL = 0;
0064       double tmp_ttrig = 0., tmp_t0 = 0., tmp_kfact = 0.;
0065       int tmp_a = 0, tmp_b = 0, tmp_c = 0, tmp_d = 0;
0066       double tmp_v_eff = 0.;
0067       while (!fileSlopes.eof()) {
0068         fileSlopes >> tmp_wheel >> tmp_sector >> tmp_station >> tmp_SL >> tmp_a >> tmp_b >> tmp_ttrig >> tmp_t0 >>
0069             tmp_kfact >> tmp_c >> tmp_d >> tmp_v_eff;
0070         vDriftEff_[tmp_wheel + 2][tmp_sector - 1][tmp_station - 1][tmp_SL - 1] = -tmp_v_eff;
0071       }
0072       fileSlopes.close();
0073     }
0074
0075     bool debug = pset.getUntrackedParameter<bool>("debug", false);
0076     fitter_ = new DTResidualFitter(debug);
0077   }
0078
0079   DTTTrigResidualCorrection::~DTTTrigResidualCorrection() {
0080     delete rootFile_;
0081     delete fitter_;
0082   }
0083
0084   void DTTTrigResidualCorrection::setES(const EventSetup& setup) {
0085     // Get tTrig record from DB
0086     ESHandle<DTTtrig> tTrig;
0087     tTrig = setup.getHandle(ttrigToken_);
0088     tTrigMap_ = &*tTrig;
0089
0090     // Get vDrift record
0091     if (readLegacyVDriftDB) {
0092       ESHandle<DTMtime> mTimeHandle;
0093       mTimeMap_ = &setup.getData(mTimeMapToken_);
0094       vDriftMap_ = nullptr;
0095     } else {
0096       ESHandle<DTRecoConditions> hVdrift;
0097       vDriftMap_ = &setup.getData(vDriftMapToken_);
0098       mTimeMap_ = nullptr;
0099       // Consistency check: no parametrization is implemented for the time being
0100       int version = vDriftMap_->version();
0101       if (version != 1) {
0102         throw cms::Exception("Configuration") << "only version 1 is presently supported for VDriftDB";
0103       }
0104     }
0105   }
0106
0107   DTTTrigData DTTTrigResidualCorrection::correction(const DTSuperLayerId& slId) {
0108     float tTrigMean, tTrigSigma, kFactor;
0109     int status = tTrigMap_->get(slId, tTrigMean, tTrigSigma, kFactor, DTTimeUnits::ns);
0110     if (status != 0)
0111       throw cms::Exception("[DTTTrigResidualCorrection]") << "Could not find tTrig entry in DB for" << slId << endl;
0112
0113     float vDrift, hitResolution = 0.;
0114     if (readLegacyVDriftDB) {  // Legacy format
0115       status = mTimeMap_->get(slId, vDrift, hitResolution, DTVelocityUnits::cm_per_ns);
0116       if (status != 0)
0117         throw cms::Exception("[DTTTrigResidualCorrection]") << "Could not find vDrift entry in DB for" << slId << endl;
0118     } else {
0119       vDrift = vDriftMap_->get(DTWireId(slId.rawId()));
0120     }
0121
0122     TH1F residualHisto = *(getHisto(slId));
0123     LogTrace("Calibration") << "[DTTTrigResidualCorrection]: \n"
0124                             << "   Mean, RMS     = " << residualHisto.GetMean() << ", " << residualHisto.GetRMS();
0125
0126     double fitMean = -1.;
0127     double fitSigma = -1.;
0128     if (useFit_) {
0129       LogTrace("Calibration") << "[DTTTrigResidualCorrection]: Fitting histogram " << residualHisto.GetName();
0130       int nSigmas = 2;
0131
0132       DTResidualFitResult fitResult = fitter_->fitResiduals(residualHisto, nSigmas);
0133       fitMean = fitResult.fitMean;
0134       fitSigma = fitResult.fitSigma;
0135
0136       LogTrace("Calibration") << "[DTTTrigResidualCorrection]: \n"
0137                               << "   Fit Mean      = " << fitMean << "\n"
0138                               << "   Fit Sigma     = " << fitSigma;
0139     }
0140     double resMean = (useFit_) ? fitMean : residualHisto.GetMean();
0141
0142     int wheel = slId.wheel();
0143     int sector = slId.sector();
0144     int station = slId.station();
0145     int superLayer = slId.superLayer();
0146     double resTime = 0.;
0147     if (useSlopesCalib_) {
0148       double vdrift_eff = vDriftEff_[wheel + 2][sector - 1][station - 1][superLayer - 1];
0149       if (vdrift_eff == 0)
0150         vdrift_eff = vDrift;
0151
0152       if (vdrift_eff)
0153         resTime = resMean / vdrift_eff;
0154
0155       LogTrace("Calibration") << "[DTTTrigResidualCorrection]: Effective vDrift, correction to tTrig = " << vdrift_eff
0156                               << ", " << resTime;
0157     } else {
0158       if (vDrift)
0159         resTime = resMean / vDrift;
0160
0161       LogTrace("Calibration") << "[DTTTrigResidualCorrection]: vDrift from DB, correction to tTrig = " << vDrift << ", "
0162                               << resTime;
0163     }
0164
0165     double corrMean = tTrigMean;
0166     double corrSigma = (tTrigSigma != 0.) ? tTrigSigma : 1.;
0167     double corrKFact = (tTrigSigma != 0.) ? (kFactor + resTime / tTrigSigma) : resTime;
0168
0169     return DTTTrigData(corrMean, corrSigma, corrKFact);
0170   }
0171
0172   const TH1F* DTTTrigResidualCorrection::getHisto(const DTSuperLayerId& slId) {
0173     string histoName = getHistoName(slId);
0174     LogTrace("Calibration") << "[DTTTrigResidualCorrection]: Accessing histogram " << histoName.c_str();
0175     TH1F* histo = static_cast<TH1F*>(rootFile_->Get(histoName.c_str()));
0176     if (!histo)
0177       throw cms::Exception("[DTTTrigResidualCorrection]") << "residual histogram not found:" << histoName << endl;
0178     return histo;
0179   }
0180
0181   string DTTTrigResidualCorrection::getHistoName(const DTSuperLayerId& slId) {
0182     int step = 3;
0183
0184     std::string wheel = std::to_string(slId.wheel());
0185     std::string station = std::to_string(slId.station());
0186     std::string sector = std::to_string(slId.sector());
0187     std::string superLayer = std::to_string(slId.superlayer());
0188     std::string Step = std::to_string(step);
0189
0190     string histoName = rootBaseDir_ + "/Wheel" + wheel + "/Station" + station + "/Sector" + sector + "/hResDist_STEP" +
0191                        Step + "_W" + wheel + "_St" + station + "_Sec" + sector + "_SL" + superLayer;
0192
0193     return histoName;
0194   }
0195
0196 }  // namespace dtCalibration