test/DBTools/DTVDriftAnalyzer.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author S. Bolognesi - INFN Torino
0005  */
0006
0007 #include "DTVDriftAnalyzer.h"
0008
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "CondFormats/DTObjects/interface/DTMtime.h"
0012 #include "CondFormats/DataRecord/interface/DTMtimeRcd.h"
0013 #include "CondFormats/DTObjects/interface/DTRecoConditions.h"
0014 #include "CondFormats/DataRecord/interface/DTRecoConditionsVdriftRcd.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include <iostream>
0017 #include "TFile.h"
0018 #include "TH1D.h"
0019 #include "TString.h"
0020
0021 using namespace edm;
0022 using namespace std;
0023
0024 DTVDriftAnalyzer::DTVDriftAnalyzer(const ParameterSet& pset)
0025     : readLegacyVDriftDB(pset.getParameter<bool>("readLegacyVDriftDB")) {
0026   // The root file which will contain the histos
0027   string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0028   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0029   theFile->cd();
0030   mTimeMapToken_ = esConsumes<edm::Transition::BeginRun>();
0031   vDriftMapToken_ = esConsumes<edm::Transition::BeginRun>();
0032 }
0033
0034 DTVDriftAnalyzer::~DTVDriftAnalyzer() { theFile->Close(); }
0035
0036 void DTVDriftAnalyzer::beginRun(const edm::Run& run, const edm::EventSetup& eventSetup) {
0037   edm::LogInfo("DTVDriftAnalyzer") << "doing beginRun" << endl;
0038   if (readLegacyVDriftDB) {
0039     ESHandle<DTMtime> mTime = eventSetup.getHandle(mTimeMapToken_);
0040     mTimeMap = &*mTime;
0041     //ESHandle<DTMtime> mTimeHandle;
0042     //mTimeMap = &eventSetup.getData(mTimeMapToken_);
0043     //vDriftMap_ = nullptr;
0044     edm::LogInfo("DTVDriftAnalyzer") << "[DTVDriftAnalyzer] MTime version: " << mTime->version() << endl;
0045   } else {
0046     ESHandle<DTRecoConditions> hVdrift = eventSetup.getHandle(vDriftMapToken_);
0047     vDriftMap_ = &*hVdrift;
0048     mTimeMap = nullptr;
0049     // Consistency check: no parametrization is implemented for the time being
0050     int version = vDriftMap_->version();
0051     if (version != 1) {
0052       throw cms::Exception("Configuration") << "only version 1 is presently supported for VDriftDB";
0053     }
0054   }
0055 }
0056
0057 void DTVDriftAnalyzer::endJob() {
0058   // Loop over DB entries
0059
0060   map<uint32_t, pair<float, float>> values;
0061
0062   if (readLegacyVDriftDB) {
0063     edm::LogInfo("DTVDriftAnalyzer") << "Reading Legacy VDrift DB" << endl;
0064
0065     for (DTMtime::const_iterator mtime = mTimeMap->begin(); mtime != mTimeMap->end(); ++mtime) {
0066       //edm::LogInfo("DTVDrift")<< typeid(mtime).name() <<endl;
0067
0068       DTWireId wireId(
0069           (*mtime).first.wheelId, (*mtime).first.stationId, (*mtime).first.sectorId, (*mtime).first.slId, 0, 0);
0070       float vdrift;
0071       float reso;
0072       DetId detId(wireId.rawId());
0073       // vdrift is cm/ns , resolution is cm
0074       mTimeMap->get(detId, vdrift, reso, DTVelocityUnits::cm_per_ns);
0075       values[wireId.rawId()] = make_pair(vdrift, reso);
0076     }
0077
0078   } else {
0079     for (DTRecoConditions::const_iterator vd = vDriftMap_->begin(); vd != vDriftMap_->end(); ++vd) {
0080       DTWireId wireId(vd->first);
0081       float vdrift = vDriftMap_->get(wireId);
0082       values[vd->first] = make_pair(vdrift, 0.f);
0083     }
0084   }
0085
0086   for (map<uint32_t, pair<float, float>>::const_iterator it = values.begin(); it != values.end(); ++it) {
0087     float vdrift = it->second.first;
0088     float reso = it->second.second;
0089     DTWireId wireId(it->first);
0090     if (wireId.wheel() == 0 && wireId.superlayer() == 1)
0091       // Print only for wheel=0 and SL=1
0092       // vdrift is cm/ns , resolution is cm
0093       edm::LogVerbatim("DTVDriftAnalyzer") << "Wire: " << wireId << endl
0094                                            << " vdrift (cm/ns): " << vdrift << endl
0095                                            << " reso (cm): " << reso << endl;
0096
0097     //Define an histo for each wheel and each superlayer type
0098     TH1D* hVDriftHisto = theVDriftHistoMap[make_pair(wireId.wheel(), wireId.superlayer())];
0099     TH1D* hResoHisto = theResoHistoMap[make_pair(wireId.wheel(), wireId.superlayer())];
0100     if (hVDriftHisto == 0) {
0101       theFile->cd();
0102       TString name = getHistoName(wireId).c_str();
0103       if (wireId.superlayer() != 2) {
0104         hVDriftHisto = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 50, 0, 50);
0105         hResoHisto = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 50, 0, 50);
0106       } else {
0107         hVDriftHisto = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 36, 0, 36);
0108         hResoHisto = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 36, 0, 36);
0109       }
0110       theVDriftHistoMap[make_pair(wireId.wheel(), wireId.superlayer())] = hVDriftHisto;
0111       theResoHistoMap[make_pair(wireId.wheel(), wireId.superlayer())] = hResoHisto;
0112     }
0113
0114     //Fill the histos and set the bin label
0115     int binNumber = wireId.sector() + 12 * (wireId.station() - 1);
0116     hVDriftHisto->SetBinContent(binNumber, vdrift);
0117     hResoHisto->SetBinContent(binNumber, reso);
0118     string labelName;
0119     stringstream theStream;
0120     if (wireId.sector() == 1)
0121       theStream << "MB" << wireId.station() << "_Sec" << wireId.sector();
0122     else
0123       theStream << "Sec" << wireId.sector();
0124     theStream >> labelName;
0125     hVDriftHisto->GetXaxis()->SetBinLabel(binNumber, labelName.c_str());
0126     hResoHisto->GetXaxis()->SetBinLabel(binNumber, labelName.c_str());
0127
0128     //Define a distribution for each wheel,station and each superlayer type
0129     vector<int> Wh_St_SL;
0130     Wh_St_SL.push_back(wireId.wheel());
0131     Wh_St_SL.push_back(wireId.station());
0132     Wh_St_SL.push_back(wireId.superlayer());
0133     TH1D* hVDriftDistrib = theVDriftDistribMap[Wh_St_SL];
0134     TH1D* hResoDistrib = theResoDistribMap[Wh_St_SL];
0135     if (hVDriftDistrib == 0) {
0136       theFile->cd();
0137       TString name = getDistribName(wireId).c_str();
0138       hVDriftDistrib = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 100, 0.00530, 0.00580);
0139       hResoDistrib = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 300, 0.015, 0.045);
0140       theVDriftDistribMap[Wh_St_SL] = hVDriftDistrib;
0141       theResoDistribMap[Wh_St_SL] = hResoDistrib;
0142     }
0143     //Fill the distributions
0144     hVDriftDistrib->Fill(vdrift);
0145     hResoDistrib->Fill(reso);
0146   }
0147
0148   //Write histos in a .root file
0149   theFile->cd();
0150   for (map<pair<int, int>, TH1D*>::const_iterator lHisto = theVDriftHistoMap.begin(); lHisto != theVDriftHistoMap.end();
0151        ++lHisto) {
0152     (*lHisto).second->GetXaxis()->LabelsOption("v");
0153     (*lHisto).second->Write();
0154   }
0155   for (map<pair<int, int>, TH1D*>::const_iterator lHisto = theResoHistoMap.begin(); lHisto != theResoHistoMap.end();
0156        ++lHisto) {
0157     (*lHisto).second->GetXaxis()->LabelsOption("v");
0158     (*lHisto).second->Write();
0159   }
0160   for (map<vector<int>, TH1D*>::const_iterator lDistrib = theVDriftDistribMap.begin();
0161        lDistrib != theVDriftDistribMap.end();
0162        ++lDistrib) {
0163     (*lDistrib).second->Write();
0164   }
0165   for (map<vector<int>, TH1D*>::const_iterator lDistrib = theResoDistribMap.begin();
0166        lDistrib != theResoDistribMap.end();
0167        ++lDistrib) {
0168     (*lDistrib).second->Write();
0169   }
0170 }
0171
0172 string DTVDriftAnalyzer::getHistoName(const DTWireId& wId) const {
0173   string histoName;
0174   stringstream theStream;
0175   theStream << "Wheel" << wId.wheel() << "_SL" << wId.superlayer();
0176   theStream >> histoName;
0177   return histoName;
0178 }
0179
0180 string DTVDriftAnalyzer::getDistribName(const DTWireId& wId) const {
0181   string histoName;
0182   stringstream theStream;
0183   theStream << "Wheel" << wId.wheel() << "_Station" << wId.station() << "_SL" << wId.superlayer();
0184   theStream >> histoName;
0185   return histoName;
0186 }