Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-17 23:31:52

0001 // -*- C++ -*-
0002 //
0003 // Package:    SiStripBaselineAnalyzer
0004 // Class:      SiStripBaselineAnalyzer
0005 //
0006 /**\class SiStripBaselineAnalyzer SiStripBaselineAnalyzer.cc Validation/SiStripAnalyzer/src/SiStripBaselineAnalyzer.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Ivan Amos Cali
0015 //         Created:  Mon Jul 28 14:10:52 CEST 2008
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022 
0023 // user include files
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 //ROOT inclusion
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 // class decleration
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       //std::cout << "bas id: " << itDSBaseline->id << " raw id: " << detId << std::endl;
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     //std::cout << "processing module N: " << actualModule_<< " detId: " << detId << " event: "<< event << std::endl;
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         //std::cout << "totADC " << totADC << std::endl;
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             //std::cout << "Found cluster in detId " << detId << " " << firststrip << " " << clus->amplitudes().size() << " -----------------------------------------------" << std::endl;
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 // ------------ method called once each job just before starting event loop  ------------
0381 void SiStripBaselineAnalyzer::beginJob() { actualModule_ = 0; }
0382 
0383 //define this as a plug-in
0384 DEFINE_FWK_MODULE(SiStripBaselineAnalyzer);