DTCalibration/plugins/DTTPAnalyzer.cc

0001 /** \class DTTPAnalyzer
0002  *
0003  *  \author A. Vilela Pereira
0004  */
0005
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0011 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0012
0013 #include <string>
0014 #include <map>
0015
0016 //class DTT0;
0017 class DTLayerId;
0018 class DTWireId;
0019 class DTGeometry;
0020 class DTTTrigBaseSync;
0021 class TFile;
0022
0023 class DTTPAnalyzer : public edm::one::EDAnalyzer<> {
0024 public:
0025   DTTPAnalyzer(const edm::ParameterSet&);
0026   ~DTTPAnalyzer() override;
0027
0028   void analyze(const edm::Event&, const edm::EventSetup&) override;
0029   void endJob() override;
0030
0031 private:
0032   std::string getHistoName(const DTLayerId&);
0033
0034   const bool subtractT0_;
0035   const edm::EDGetTokenT<DTDigiCollection> digiToken_;
0036
0037   TFile* rootFile_;
0038   edm::ESHandle<DTGeometry> dtGeom_;
0039   const edm::ESGetToken<DTGeometry, MuonGeometryRecord> dtGeomToken_;
0040   std::unique_ptr<DTTTrigBaseSync> tTrigSync_;
0041
0042   // Map of the t0 and sigma histos by layer
0043   std::map<DTWireId, int> nDigisPerWire_;
0044   std::map<DTWireId, double> sumWPerWire_;
0045   std::map<DTWireId, double> sumW2PerWire_;
0046 };
0047
0048 #include "FWCore/Framework/interface/Event.h"
0049 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0050 #include "FWCore/Framework/interface/MakerMacros.h"
0051
0052 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
0053 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0054
0055 #include "Geometry/Records/interface/MuonNumberingRecord.h"
0056 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0057
0058 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0059
0060 #include "CondFormats/DTObjects/interface/DTT0.h"
0061
0062 #include "TH1F.h"
0063 #include "TFile.h"
0064
0065 DTTPAnalyzer::DTTPAnalyzer(const edm::ParameterSet& pset)
0066     : subtractT0_(pset.getParameter<bool>("subtractT0")),
0067       digiToken_(consumes<DTDigiCollection>(pset.getParameter<edm::InputTag>("digiLabel"))),
0068       dtGeomToken_(esConsumes()) {
0069   std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName");
0070   rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
0071   rootFile_->cd();
0072
0073   if (subtractT0_)
0074     tTrigSync_ = DTTTrigSyncFactory::get()->create(pset.getParameter<std::string>("tTrigMode"),
0075                                                    pset.getParameter<edm::ParameterSet>("tTrigModeConfig"),
0076                                                    consumesCollector());
0077 }
0078
0079 DTTPAnalyzer::~DTTPAnalyzer() { rootFile_->Close(); }
0080
0081 void DTTPAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0082   // Get the digis from the event
0083   const edm::Handle<DTDigiCollection>& digis = event.getHandle(digiToken_);
0084
0085   if (subtractT0_) {
0086     tTrigSync_->setES(setup);
0087   }
0088   // Get the DT Geometry
0089   dtGeom_ = setup.getHandle(dtGeomToken_);
0090
0091   // Iterate through all digi collections ordered by LayerId
0092   DTDigiCollection::DigiRangeIterator dtLayerIt;
0093   for (dtLayerIt = digis->begin(); dtLayerIt != digis->end(); ++dtLayerIt) {
0094     // Get the iterators over the digis associated with this LayerId
0095     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
0096
0097     // Get the layerId
0098     const DTLayerId layerId = (*dtLayerIt).first;  //FIXME: check to be in the right sector
0099
0100     // Loop over all digis in the given layer
0101     for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; ++digi) {
0102       const DTWireId wireId(layerId, (*digi).wire());
0103
0104       double t0 = (*digi).countsTDC();
0105
0106       //FIXME: Reject digis not coming from TP
0107
0108       if (subtractT0_) {
0109         const DTLayer* layer = nullptr;  //fake
0110         const GlobalPoint glPt;          //fake
0111         double offset = tTrigSync_->offset(layer, wireId, glPt);
0112         t0 -= offset;
0113       }
0114
0115       if (nDigisPerWire_.find(wireId) == nDigisPerWire_.end()) {
0116         nDigisPerWire_[wireId] = 0;
0117         sumWPerWire_[wireId] = 0.;
0118         sumW2PerWire_[wireId] = 0.;
0119       }
0120
0121       ++nDigisPerWire_[wireId];
0122       sumWPerWire_[wireId] += t0;
0123       sumW2PerWire_[wireId] += t0 * t0;
0124     }
0125   }
0126 }
0127
0128 void DTTPAnalyzer::endJob() {
0129   rootFile_->cd();
0130   std::map<DTLayerId, TH1F*> meanHistoMap;
0131   std::map<DTLayerId, TH1F*> sigmaHistoMap;
0132   for (std::map<DTWireId, int>::const_iterator wireIdIt = nDigisPerWire_.begin(); wireIdIt != nDigisPerWire_.end();
0133        ++wireIdIt) {
0134     DTWireId wireId((*wireIdIt).first);
0135
0136     int nDigis = nDigisPerWire_[wireId];
0137     double sumW = sumWPerWire_[wireId];
0138     double sumW2 = sumW2PerWire_[wireId];
0139
0140     double mean = sumW / nDigis;
0141     double rms = sumW2 / nDigis - mean * mean;
0142     rms = sqrt(rms);
0143
0144     DTLayerId layerId = wireId.layerId();
0145     if (meanHistoMap.find(layerId) == meanHistoMap.end()) {
0146       std::string histoName = getHistoName(layerId);
0147       const int firstChannel = dtGeom_->layer(layerId)->specificTopology().firstChannel();
0148       const int nWires = dtGeom_->layer(layerId)->specificTopology().channels();
0149       TH1F* meanHistoTP = new TH1F((histoName + "_tpMean").c_str(),
0150                                    "mean from test pulses by channel",
0151                                    nWires,
0152                                    firstChannel,
0153                                    (firstChannel + nWires));
0154       TH1F* sigmaHistoTP = new TH1F((histoName + "_tpSigma").c_str(),
0155                                     "sigma from test pulses by channel",
0156                                     nWires,
0157                                     firstChannel,
0158                                     (firstChannel + nWires));
0159       meanHistoMap[layerId] = meanHistoTP;
0160       sigmaHistoMap[layerId] = sigmaHistoTP;
0161     }
0162     // Fill the histograms
0163     int nBin = meanHistoMap[layerId]->GetXaxis()->FindFixBin(wireId.wire());
0164     meanHistoMap[layerId]->SetBinContent(nBin, mean);
0165     sigmaHistoMap[layerId]->SetBinContent(nBin, rms);
0166   }
0167
0168   for (std::map<DTLayerId, TH1F*>::const_iterator key = meanHistoMap.begin(); key != meanHistoMap.end(); ++key) {
0169     meanHistoMap[(*key).first]->Write();
0170     sigmaHistoMap[(*key).first]->Write();
0171   }
0172 }
0173
0174 std::string DTTPAnalyzer::getHistoName(const DTLayerId& lId) {
0175   std::string histoName;
0176   std::stringstream theStream;
0177   theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector() << "_SL" << lId.superlayer() << "_L"
0178             << lId.layer();
0179   theStream >> histoName;
0180   return histoName;
0181 }
0182
0183 DEFINE_FWK_MODULE(DTTPAnalyzer);