File indexing completed on 2024-04-06 12:26:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <iostream>
0022
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "DataFormats/Common/interface/DetSet.h"
0032 #include "DataFormats/Common/interface/DetSetVector.h"
0033 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0034
0035 #include "DataFormats/SiStripDigi/interface/SiStripProcessedRawDigi.h"
0036 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
0037 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0038 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
0039
0040 #include "CondFormats/SiStripObjects/interface/SiStripPedestals.h"
0041 #include "CondFormats/DataRecord/interface/SiStripPedestalsRcd.h"
0042
0043 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripPedestalsSubtractor.h"
0044 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripCommonModeNoiseSubtractor.h"
0045 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingFactory.h"
0046
0047 #include "FWCore/ServiceRegistry/interface/Service.h"
0048 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0049 #include "CommonTools/Utils/interface/TFileDirectory.h"
0050
0051 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0052
0053
0054 #include "TROOT.h"
0055 #include "TFile.h"
0056 #include "TNtuple.h"
0057 #include "TMath.h"
0058 #include "TCanvas.h"
0059 #include "TH1F.h"
0060 #include "TH2F.h"
0061 #include "TProfile.h"
0062 #include "TList.h"
0063 #include "TString.h"
0064 #include "TStyle.h"
0065 #include "TGraph.h"
0066 #include "TMultiGraph.h"
0067 #include "THStack.h"
0068
0069
0070
0071
0072
0073 class SiStripBaselineAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0074 public:
0075 explicit SiStripBaselineAnalyzer(const edm::ParameterSet&);
0076 ~SiStripBaselineAnalyzer() override;
0077
0078 private:
0079 void beginJob() override;
0080 void analyze(const edm::Event&, const edm::EventSetup&) override;
0081
0082 std::unique_ptr<SiStripPedestalsSubtractor> subtractorPed_;
0083 edm::ESGetToken<SiStripPedestals, SiStripPedestalsRcd> pedestalsToken_;
0084 const SiStripPedestals* pedestalsHandle;
0085 edm::ESWatcher<SiStripPedestalsRcd> pedestalsWatcher_;
0086 std::vector<int> pedestals;
0087
0088 bool plotClusters_;
0089 bool plotBaseline_;
0090 bool plotBaselinePoints_;
0091 bool plotRawDigi_;
0092 bool plotDigis_;
0093 bool plotAPVCM_;
0094 bool plotPedestals_;
0095
0096 edm::InputTag srcBaseline_;
0097 edm::InputTag srcBaselinePoints_;
0098 edm::InputTag srcAPVCM_;
0099 edm::InputTag srcProcessedRawDigi_;
0100 edm::InputTag srcDigis_;
0101
0102 edm::EDGetTokenT<edm::DetSetVector<SiStripProcessedRawDigi>> APVCMToken;
0103 edm::EDGetTokenT<edm::DetSetVector<SiStripRawDigi>> processedRawToken;
0104 edm::EDGetTokenT<edm::DetSetVector<SiStripProcessedRawDigi>> baselineToken;
0105 edm::EDGetTokenT<edm::DetSetVector<SiStripDigi>> baselinePointsToken;
0106 edm::EDGetTokenT<edm::DetSetVector<SiStripDigi>> digisToken;
0107 edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> clustersToken;
0108
0109 edm::Service<TFileService> fs_;
0110
0111 TH1F* h1BadAPVperEvent_;
0112
0113 TH1F* h1ProcessedRawDigis_;
0114 TH1F* h1Baseline_;
0115 TH1F* h1Clusters_;
0116 TH1F* h1APVCM_;
0117 TH1F* h1Pedestals_;
0118
0119 TCanvas* Canvas_;
0120 std::vector<TH1F> vProcessedRawDigiHisto_;
0121 std::vector<TH1F> vBaselineHisto_;
0122 std::vector<TH1F> vBaselinePointsHisto_;
0123 std::vector<TH1F> vClusterHisto_;
0124
0125 uint16_t nModuletoDisplay_;
0126 uint16_t actualModule_;
0127 };
0128
0129 SiStripBaselineAnalyzer::SiStripBaselineAnalyzer(const edm::ParameterSet& conf) {
0130 usesResource("TFileService");
0131
0132 pedestalsToken_ = esConsumes();
0133 srcBaseline_ = conf.getParameter<edm::InputTag>("srcBaseline");
0134 srcBaselinePoints_ = conf.getParameter<edm::InputTag>("srcBaselinePoints");
0135 srcProcessedRawDigi_ = conf.getParameter<edm::InputTag>("srcProcessedRawDigi");
0136 srcDigis_ = conf.getParameter<edm::InputTag>("srcDigis");
0137 srcAPVCM_ = conf.getParameter<edm::InputTag>("srcAPVCM");
0138
0139 APVCMToken = consumes<edm::DetSetVector<SiStripProcessedRawDigi>>(srcAPVCM_);
0140 processedRawToken = consumes<edm::DetSetVector<SiStripRawDigi>>(srcProcessedRawDigi_);
0141 baselineToken = consumes<edm::DetSetVector<SiStripProcessedRawDigi>>(srcBaseline_);
0142 baselinePointsToken = consumes<edm::DetSetVector<SiStripDigi>>(srcBaselinePoints_);
0143 digisToken = consumes<edm::DetSetVector<SiStripDigi>>(srcDigis_);
0144 edm::InputTag clusLabel("siStripClusters");
0145 clustersToken = consumes<edmNew::DetSetVector<SiStripCluster>>(clusLabel);
0146
0147 subtractorPed_ = SiStripRawProcessingFactory::create_SubtractorPed(conf.getParameter<edm::ParameterSet>("Algorithms"),
0148 consumesCollector());
0149 nModuletoDisplay_ = conf.getParameter<uint32_t>("nModuletoDisplay");
0150 plotClusters_ = conf.getParameter<bool>("plotClusters");
0151 plotBaseline_ = conf.getParameter<bool>("plotBaseline");
0152 plotBaselinePoints_ = conf.getParameter<bool>("plotBaselinePoints");
0153 plotRawDigi_ = conf.getParameter<bool>("plotRawDigi");
0154 plotAPVCM_ = conf.getParameter<bool>("plotAPVCM");
0155 plotPedestals_ = conf.getParameter<bool>("plotPedestals");
0156 plotDigis_ = conf.getParameter<bool>("plotDigis");
0157
0158 h1BadAPVperEvent_ = fs_->make<TH1F>("BadAPV/Event", "BadAPV/Event", 2001, -0.5, 2000.5);
0159 h1BadAPVperEvent_->SetXTitle("# Modules with Bad APVs");
0160 h1BadAPVperEvent_->SetYTitle("Entries");
0161 h1BadAPVperEvent_->SetLineWidth(2);
0162 h1BadAPVperEvent_->SetLineStyle(2);
0163
0164 h1APVCM_ = fs_->make<TH1F>("APV CM", "APV CM", 2048, -1023.5, 1023.5);
0165 h1APVCM_->SetXTitle("APV CM [adc]");
0166 h1APVCM_->SetYTitle("Entries");
0167 h1APVCM_->SetLineWidth(2);
0168 h1APVCM_->SetLineStyle(2);
0169
0170 h1Pedestals_ = fs_->make<TH1F>("Pedestals", "Pedestals", 2048, -1023.5, 1023.5);
0171 h1Pedestals_->SetXTitle("Pedestals [adc]");
0172 h1Pedestals_->SetYTitle("Entries");
0173 h1Pedestals_->SetLineWidth(2);
0174 h1Pedestals_->SetLineStyle(2);
0175 }
0176
0177 SiStripBaselineAnalyzer::~SiStripBaselineAnalyzer() = default;
0178
0179 void SiStripBaselineAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& es) {
0180 using namespace edm;
0181 if (plotPedestals_ && actualModule_ == 0) {
0182 if (pedestalsWatcher_.check(es)) {
0183 pedestalsHandle = &es.getData(pedestalsToken_);
0184 }
0185
0186 std::vector<uint32_t> detIdV;
0187 pedestalsHandle->getDetIds(detIdV);
0188
0189 for (uint32_t i = 0; i < detIdV.size(); ++i) {
0190 pedestals.clear();
0191 SiStripPedestals::Range pedestalsRange = pedestalsHandle->getRange(detIdV[i]);
0192 pedestals.resize((pedestalsRange.second - pedestalsRange.first) * 8 / 10);
0193 pedestalsHandle->allPeds(pedestals, pedestalsRange);
0194 for (uint32_t it = 0; it < pedestals.size(); ++it)
0195 h1Pedestals_->Fill(pedestals[it]);
0196 }
0197 }
0198
0199 if (plotAPVCM_) {
0200 edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi>> moduleCM;
0201 edm::InputTag CMLabel("siStripZeroSuppression:APVCM");
0202 e.getByToken(APVCMToken, moduleCM);
0203
0204 edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itCMDetSetV = moduleCM->begin();
0205 for (; itCMDetSetV != moduleCM->end(); ++itCMDetSetV) {
0206 edm::DetSet<SiStripProcessedRawDigi>::const_iterator itCM = itCMDetSetV->begin();
0207 for (; itCM != itCMDetSetV->end(); ++itCM)
0208 h1APVCM_->Fill(itCM->adc());
0209 }
0210 }
0211
0212 subtractorPed_->init(es);
0213
0214 edm::Handle<edm::DetSetVector<SiStripRawDigi>> moduleRawDigi;
0215 if (plotRawDigi_)
0216 e.getByToken(processedRawToken, moduleRawDigi);
0217
0218 edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi>> moduleBaseline;
0219 if (plotBaseline_)
0220 e.getByToken(baselineToken, moduleBaseline);
0221
0222 edm::Handle<edm::DetSetVector<SiStripDigi>> moduleBaselinePoints;
0223 if (plotBaselinePoints_)
0224 e.getByToken(baselinePointsToken, moduleBaselinePoints);
0225
0226 edm::Handle<edm::DetSetVector<SiStripDigi>> moduleDigis;
0227 if (plotDigis_)
0228 e.getByToken(digisToken, moduleDigis);
0229
0230 edm::Handle<edmNew::DetSetVector<SiStripCluster>> clusters;
0231 if (plotClusters_) {
0232 e.getByToken(clustersToken, clusters);
0233 }
0234
0235 char detIds[20];
0236 char evs[20];
0237 char runs[20];
0238
0239 TFileDirectory sdProcessedRawDigis_ = fs_->mkdir("ProcessedRawDigis");
0240 TFileDirectory sdBaseline_ = fs_->mkdir("Baseline");
0241 TFileDirectory sdBaselinePoints_ = fs_->mkdir("BaselinePoints");
0242 TFileDirectory sdClusters_ = fs_->mkdir("Clusters");
0243 TFileDirectory sdDigis_ = fs_->mkdir("Digis");
0244
0245 edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itDSBaseline;
0246 if (plotBaseline_)
0247 itDSBaseline = moduleBaseline->begin();
0248 edm::DetSetVector<SiStripRawDigi>::const_iterator itRawDigis = moduleRawDigi->begin();
0249
0250 uint32_t NBabAPVs = moduleRawDigi->size();
0251 std::cout << "Number of module with HIP in this event: " << NBabAPVs << std::endl;
0252 h1BadAPVperEvent_->Fill(NBabAPVs);
0253
0254 for (; itRawDigis != moduleRawDigi->end(); ++itRawDigis) {
0255 if (actualModule_ > nModuletoDisplay_)
0256 return;
0257 uint32_t detId = itRawDigis->id;
0258
0259 if (plotBaseline_) {
0260
0261 if (itDSBaseline->id != detId) {
0262 std::cout << "Collections out of Synch. Something of fishy is going on ;-)" << std::endl;
0263 return;
0264 }
0265 }
0266
0267 actualModule_++;
0268 edm::RunNumber_t const run = e.id().run();
0269 edm::EventNumber_t const event = e.id().event();
0270
0271
0272 edm::DetSet<SiStripRawDigi>::const_iterator itRaw = itRawDigis->begin();
0273 bool restAPV[6] = {false, false, false, false, false, false};
0274 int strip = 0, totADC = 0;
0275 int minAPVRes = 7, maxAPVRes = -1;
0276 for (; itRaw != itRawDigis->end(); ++itRaw, ++strip) {
0277 float adc = itRaw->adc();
0278 totADC += adc;
0279 if (strip % 127 == 0) {
0280
0281 int APV = strip / 128;
0282 if (totADC != 0) {
0283 restAPV[APV] = true;
0284 totADC = 0;
0285 if (APV > maxAPVRes)
0286 maxAPVRes = APV;
0287 if (APV < minAPVRes)
0288 minAPVRes = APV;
0289 }
0290 }
0291 }
0292
0293 uint16_t bins = 768;
0294 float minx = -0.5, maxx = 767.5;
0295 if (minAPVRes != 7) {
0296 minx = minAPVRes * 128 - 0.5;
0297 maxx = maxAPVRes * 128 + 127.5;
0298 bins = maxx - minx;
0299 }
0300
0301 sprintf(detIds, "%ul", detId);
0302 sprintf(evs, "%llu", event);
0303 sprintf(runs, "%u", run);
0304 char* dHistoName = Form("Id:%s_run:%s_ev:%s", detIds, runs, evs);
0305 h1ProcessedRawDigis_ = sdProcessedRawDigis_.make<TH1F>(dHistoName, dHistoName, bins, minx, maxx);
0306
0307 if (plotBaseline_) {
0308 h1Baseline_ = sdBaseline_.make<TH1F>(dHistoName, dHistoName, bins, minx, maxx);
0309 h1Baseline_->SetXTitle("strip#");
0310 h1Baseline_->SetYTitle("ADC");
0311 h1Baseline_->SetMaximum(1024.);
0312 h1Baseline_->SetMinimum(-300.);
0313 h1Baseline_->SetLineWidth(2);
0314 h1Baseline_->SetLineStyle(2);
0315 h1Baseline_->SetLineColor(2);
0316 }
0317
0318 if (plotClusters_) {
0319 h1Clusters_ = sdClusters_.make<TH1F>(dHistoName, dHistoName, bins, minx, maxx);
0320
0321 h1Clusters_->SetXTitle("strip#");
0322 h1Clusters_->SetYTitle("ADC");
0323 h1Clusters_->SetMaximum(1024.);
0324 h1Clusters_->SetMinimum(-300.);
0325 h1Clusters_->SetLineWidth(2);
0326 h1Clusters_->SetLineStyle(2);
0327 h1Clusters_->SetLineColor(3);
0328 }
0329
0330 h1ProcessedRawDigis_->SetXTitle("strip#");
0331 h1ProcessedRawDigis_->SetYTitle("ADC");
0332 h1ProcessedRawDigis_->SetMaximum(1024.);
0333 h1ProcessedRawDigis_->SetMinimum(-300.);
0334 h1ProcessedRawDigis_->SetLineWidth(2);
0335
0336 std::vector<int16_t> ProcessedRawDigis(itRawDigis->size());
0337 subtractorPed_->subtract(*itRawDigis, ProcessedRawDigis);
0338
0339 edm::DetSet<SiStripProcessedRawDigi>::const_iterator itBaseline;
0340 if (plotBaseline_)
0341 itBaseline = itDSBaseline->begin();
0342 std::vector<int16_t>::const_iterator itProcessedRawDigis;
0343
0344 strip = 0;
0345 for (itProcessedRawDigis = ProcessedRawDigis.begin(); itProcessedRawDigis != ProcessedRawDigis.end();
0346 ++itProcessedRawDigis) {
0347 if (restAPV[strip / 128]) {
0348 float adc = *itProcessedRawDigis;
0349 h1ProcessedRawDigis_->Fill(strip, adc);
0350 if (plotBaseline_) {
0351 h1Baseline_->Fill(strip, itBaseline->adc());
0352 ++itBaseline;
0353 }
0354 }
0355 ++strip;
0356 }
0357
0358 if (plotBaseline_)
0359 ++itDSBaseline;
0360 if (plotClusters_) {
0361 edmNew::DetSetVector<SiStripCluster>::const_iterator itClusters = clusters->begin();
0362 for (; itClusters != clusters->end(); ++itClusters) {
0363 for (edmNew::DetSet<SiStripCluster>::const_iterator clus = itClusters->begin(); clus != itClusters->end();
0364 ++clus) {
0365 if (itClusters->id() == detId) {
0366 int firststrip = clus->firstStrip();
0367
0368 strip = 0;
0369 for (auto itAmpl = clus->amplitudes().begin(); itAmpl != clus->amplitudes().end(); ++itAmpl) {
0370 h1Clusters_->Fill(firststrip + strip, *itAmpl);
0371 ++strip;
0372 }
0373 }
0374 }
0375 }
0376 }
0377 }
0378 }
0379
0380
0381 void SiStripBaselineAnalyzer::beginJob() { actualModule_ = 0; }
0382
0383
0384 DEFINE_FWK_MODULE(SiStripBaselineAnalyzer);