SiStrip/plugins/SiStripApvGainInspector.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CondTools/SiStrip
0004 // Class:      SiStripApvGainInspector
0005 //
0006 /*
0007  *\class SiStripApvGainInspector SiStripApvGainInspector.cc CalibTracker/SiStripChannelGain/plugins/SiStripApvGainInspector.cc
0008
0009  Description: This module allows redo the per-APV gain fits with different PDFs (landau, landau + gaus convolution, etc.) starting from the Charge vs APV index plot produced in the SiStrip G2 APV gain PCL workflow. It is possible to inspect the 1D charge distributions for certain APVs after fitting by means of specifying them via the parameter selectedModules.
0010
0011  Implementation: largely based off CalibTracker/SiStripChannelGain/src/SiStripGainsPCLHarvester.cc
0012
0013 */
0014 //
0015 // Original Author:  Marco Musich
0016 //         Created:  Tue, 05 Jun 2018 15:46:15 GMT
0017 //
0018 //
0019
0020 // system include files
0021 #include <cmath> /* log */
0022 #include <memory>
0023
0024 // user include files
0025 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0026 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0027 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0028 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
0029 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
0030 #include "CalibTracker/SiStripChannelGain/interface/APVGainStruct.h"
0031 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0032 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0035 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
0036 #include "CondTools/SiStrip/interface/SiStripMiscalibrateHelper.h"
0037 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0038 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0039 #include "FWCore/Framework/interface/ESHandle.h"
0040 #include "FWCore/Framework/interface/Event.h"
0041 #include "FWCore/Framework/interface/EventSetup.h"
0042 #include "FWCore/Framework/interface/Frameworkfwd.h"
0043 #include "FWCore/Framework/interface/MakerMacros.h"
0044 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0045 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "FWCore/Utilities/interface/InputTag.h"
0048 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0049 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0051 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0052 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0053 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0054 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0055
0056 // ROOT includes
0057 #include "TStyle.h"
0058 #include "TCanvas.h"
0059 #include "TFile.h"
0060 #include "TTree.h"
0061 #include "TH1F.h"
0062 #include "TH2S.h"
0063 #include "TProfile.h"
0064 #include "TF1.h"
0065
0066 //
0067 // class declaration
0068 //
0069 class SiStripApvGainInspector : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0070 public:
0071   explicit SiStripApvGainInspector(const edm::ParameterSet&);
0072   ~SiStripApvGainInspector() override;
0073
0074   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0075
0076 private:
0077   void beginJob() override;
0078   void analyze(const edm::Event&, const edm::EventSetup&) override;
0079   void endJob() override;
0080   void checkBookAPVColls(const edm::EventSetup& es);
0081   void checkAndRetrieveTopology(const edm::EventSetup& setup);
0082   bool isGoodLandauFit(double* FitResults);
0083   void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange = 50, double HighRange = 5400);
0084   void getPeakOfLanGau(TH1* InputHisto, double* FitResults, double LowRange = 50, double HighRange = 5400);
0085   void doFakeFit(TH1* InputHisto, double* FitResults);
0086   void getPeakOfLandauAroundMax(TH1* InputHisto, double* FitResults, double LowRange = 100, double HighRange = 100);
0087   static double langaufun(Double_t* x, Double_t* par);
0088   void storeOnTree(TFileService* tfs);
0089   void makeNicePlotStyle(TH1F* plot);
0090   std::unique_ptr<SiStripApvGain> getNewObject();
0091   std::map<std::string, TH1*> bookQualityMonitor(const TFileDirectory& dir);
0092   void fillQualityMonitor();
0093
0094   void inline fill1D(std::map<std::string, TH1*>& h, const std::string& s, double x) {
0095     if (h.count(s) == 0) {
0096       edm::LogWarning("SiStripApvGainInspector") << "Trying to fill non-existing Histogram named " << s << std::endl;
0097       return;
0098     }
0099     h[s]->Fill(x);
0100   }
0101
0102   void inline fill2D(std::map<std::string, TH1*>& h, const std::string& s, double x, double y) {
0103     if (h.count(s) == 0) {
0104       edm::LogWarning("SiStripApvGainInspector") << "Trying to fill non-existing Histogram named " << s << std::endl;
0105       return;
0106     }
0107     h[s]->Fill(x, y);
0108   }
0109
0110   // ----------member data ---------------------------
0111   enum fitMode { landau = 1, landauAroundMax = 2, landauGauss = 3, fake = 4 };
0112
0113   const std::vector<std::string> fitModeStrings = {
0114       "",  // Enum values start from 1, so index 0 is empty or can be used as "invalid"
0115       "landau",
0116       "landauAroundMax",
0117       "landauGauss",
0118       "fake"};
0119
0120   inline bool isValidMode(int mode) const {
0121     return mode == landau || mode == landauAroundMax || mode == landauGauss || mode == fake;
0122   }
0123
0124   const edm::ESGetToken<SiStripGain, SiStripGainRcd> gainToken_;
0125   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> qualityToken_;
0126   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0127   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0128
0129   TFileService* tfs;
0130
0131   // map the APV ids to the charge plots
0132   std::map<std::pair<unsigned char, uint32_t>, TH1F*> histoMap_;
0133
0134   edm::ESHandle<TrackerGeometry> tkGeom_;
0135   const TrackerGeometry* bareTkGeomPtr_;  // ugly hack to fill APV colls only once, but checks
0136   const TrackerTopology* tTopo_;
0137
0138   int NStripAPVs;
0139   int NPixelDets;
0140
0141   unsigned int GOOD;
0142   unsigned int BAD;
0143   unsigned int MASKED;
0144
0145   std::vector<std::shared_ptr<stAPVGain>> APVsCollOrdered;
0146   std::unordered_map<unsigned int, std::shared_ptr<stAPVGain>> APVsColl;
0147
0148   const TH2F* Charge_Vs_Index;
0149   TFile* fin;
0150   fitMode fitMode_;  // Declare the enum variable
0151   const std::string filename_;
0152   double minNrEntries;
0153   std::vector<unsigned int> wantedmods;
0154
0155   std::unique_ptr<TrackerMap> ratio_map;
0156   std::unique_ptr<TrackerMap> old_payload_map;
0157   std::unique_ptr<TrackerMap> new_payload_map;
0158   std::unique_ptr<TrackerMap> mpv_map;
0159   std::unique_ptr<TrackerMap> mpv_err_map;
0160   std::unique_ptr<TrackerMap> entries_map;
0161   std::unique_ptr<TrackerMap> fitChi2_map;
0162
0163   std::map<std::string, TH1*> hControl;
0164 };
0165
0166 //
0167 // constructors and destructor
0168 //
0169 SiStripApvGainInspector::SiStripApvGainInspector(const edm::ParameterSet& iConfig)
0170     : gainToken_(esConsumes()),
0171       qualityToken_(esConsumes()),
0172       tkGeomToken_(esConsumes()),
0173       tTopoToken_(esConsumes()),
0174       bareTkGeomPtr_(nullptr),
0175       tTopo_(nullptr),
0176       GOOD(0),
0177       BAD(0),
0178       filename_(iConfig.getUntrackedParameter<std::string>("inputFile")),
0179       minNrEntries(iConfig.getUntrackedParameter<double>("minNrEntries", 20)),
0180       wantedmods(iConfig.getUntrackedParameter<std::vector<unsigned int>>("selectedModules")) {
0181   usesResource(TFileService::kSharedResource);
0182   usesResource(cond::service::PoolDBOutputService::kSharedResource);
0183
0184   sort(wantedmods.begin(), wantedmods.end());
0185
0186   edm::LogInfo("SelectedModules") << "Selected module list";
0187   for (std::vector<unsigned int>::const_iterator mod = wantedmods.begin(); mod != wantedmods.end(); mod++) {
0188     edm::LogVerbatim("SelectedModules") << *mod;
0189   }
0190
0191   int modeValue = iConfig.getParameter<int>("fitMode");
0192   if (!isValidMode(modeValue)) {
0193     throw std::invalid_argument("Invalid value provided for 'fitMode'");
0194   } else {
0195     edm::LogPrint("SiStripApvGainInspector") << "Chosen fitting mode: " << fitModeStrings[modeValue];
0196   }
0197
0198   fitMode_ = static_cast<fitMode>(modeValue);
0199
0200   //now do what ever initialization is needed
0201   fin = TFile::Open(filename_.c_str(), "READ");
0202   Charge_Vs_Index = (TH2F*)fin->Get("DQMData/Run 999999/AlCaReco/Run summary/SiStripGainsAAG/Charge_Vs_Index_AagBunch");
0203
0204   ratio_map = std::make_unique<TrackerMap>("ratio");
0205   ratio_map->setTitle("Average by module of the G2 Gain payload ratio (new/old)");
0206   ratio_map->setPalette(1);
0207
0208   new_payload_map = std::make_unique<TrackerMap>("new_payload");
0209   new_payload_map->setTitle("Tracker Map of Updated G2 Gain payload averaged by module");
0210   new_payload_map->setPalette(1);
0211
0212   old_payload_map = std::make_unique<TrackerMap>("old_payload");
0213   old_payload_map->setTitle("Tracker Map of Starting G2 Gain Payload averaged by module");
0214   old_payload_map->setPalette(1);
0215
0216   // fit quality maps
0217
0218   mpv_map = std::make_unique<TrackerMap>("MPV");
0219   mpv_map->setTitle("Landau Fit MPV average value per module [ADC counts/mm]");
0220   mpv_map->setPalette(1);
0221
0222   mpv_err_map = std::make_unique<TrackerMap>("MPVerr");
0223   mpv_err_map->setTitle("Landau Fit MPV average error per module [ADC counts/mm]");
0224   mpv_err_map->setPalette(1);
0225
0226   entries_map = std::make_unique<TrackerMap>("Entries");
0227   entries_map->setTitle("log_{10}(entries) average per module");
0228   entries_map->setPalette(1);
0229
0230   fitChi2_map = std::make_unique<TrackerMap>("FitChi2");
0231   fitChi2_map->setTitle("log_{10}(Fit #chi^{2}/ndf) average per module");
0232   fitChi2_map->setPalette(1);
0233 }
0234
0235 // do anything here that needs to be done at desctruction time
0236 // (e.g. close files, deallocate resources etc.)
0237 SiStripApvGainInspector::~SiStripApvGainInspector() {
0238   fin->Close();
0239   delete fin;
0240 }
0241
0242 //
0243 // member functions
0244 //
0245
0246 // ------------ method called for each event  ------------
0247 void SiStripApvGainInspector::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0248   using namespace edm;
0249   this->checkBookAPVColls(iSetup);  // check whether APV colls are booked and do so if not yet done
0250   this->checkAndRetrieveTopology(iSetup);
0251
0252   edm::ESHandle<SiStripGain> gainHandle = iSetup.getHandle(gainToken_);
0253   if (!gainHandle.isValid()) {
0254     edm::LogError("SiStripApvGainInspector") << "gainHandle is not valid\n";
0255     exit(0);
0256   }
0257
0258   edm::ESHandle<SiStripQuality> SiStripQuality_ = iSetup.getHandle(qualityToken_);
0259
0260   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0261     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0262
0263     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0264       continue;
0265
0266     APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId, APV->APVId);
0267
0268     if (gainHandle->getNumberOfTags() != 2) {
0269       edm::LogError("SiStripApvGainInspector") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
0270       fflush(stdout);
0271       exit(0);
0272     };
0273     float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
0274     if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
0275       edm::LogWarning("SiStripApvGainInspector") << "WARNING: ParticleGain in the global tag changed\n";
0276     APV->PreviousGain = newPreviousGain;
0277
0278     float newPreviousGainTick = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
0279     if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
0280       edm::LogWarning("SiStripApvGainInspector")
0281           << "WARNING: TickMarkGain in the global tag changed\n"
0282           << std::endl
0283           << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
0284           << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
0285           << std::endl;
0286     }
0287     APV->PreviousGainTick = newPreviousGainTick;
0288   }
0289
0290   unsigned int I = 0;
0291   TH1F* Proj = nullptr;
0292   double FitResults[6];
0293   double MPVmean = 300;
0294
0295   if (Charge_Vs_Index == nullptr) {
0296     edm::LogError("SiStripGainsPCLHarvester") << "Harvesting: could not find input histogram " << std::endl;
0297     return;
0298   }
0299
0300   printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0301   printf("Fitting Charge Distribution  :");
0302   int TreeStep = APVsColl.size() / 50;
0303
0304   for (auto it = APVsColl.begin(); it != APVsColl.end(); it++, I++) {
0305     if (I % TreeStep == 0) {
0306       printf(".");
0307       fflush(stdout);
0308     }
0309     std::shared_ptr<stAPVGain> APV = it->second;
0310     if (APV->Bin < 0)
0311       APV->Bin = Charge_Vs_Index->GetXaxis()->FindBin(APV->Index);
0312
0313     Proj = (TH1F*)(Charge_Vs_Index->ProjectionY(
0314         "", Charge_Vs_Index->GetXaxis()->FindBin(APV->Index), Charge_Vs_Index->GetXaxis()->FindBin(APV->Index), "e"));
0315     if (!Proj)
0316       continue;
0317
0318     switch (fitMode_) {
0319       case landau:
0320         getPeakOfLandau(Proj, FitResults);
0321         break;
0322       case landauAroundMax:
0323         getPeakOfLandauAroundMax(Proj, FitResults);
0324         break;
0325       case landauGauss:
0326         getPeakOfLanGau(Proj, FitResults);
0327         break;
0328       case fake:
0329         doFakeFit(Proj, FitResults);
0330         break;
0331       default:
0332         throw std::invalid_argument("Invalid value provided for 'fitMode'");
0333     }
0334
0335     APV->FitMPV = FitResults[0];
0336     APV->FitMPVErr = FitResults[1];
0337     APV->FitWidth = FitResults[2];
0338     APV->FitWidthErr = FitResults[3];
0339     APV->FitChi2 = FitResults[4];
0340     APV->FitNorm = FitResults[5];
0341     APV->NEntries = Proj->GetEntries();
0342
0343     for (const auto& mod : wantedmods) {
0344       if (mod == APV->DetId) {
0345         edm::LogInfo("ModuleFound") << " module " << mod << " found! Storing... " << std::endl;
0346         histoMap_[std::make_pair(APV->APVId, APV->DetId)] = (TH1F*)Proj->Clone(Form("hClone_%s", Proj->GetName()));
0347       }
0348     }
0349
0350     if (isGoodLandauFit(FitResults)) {
0351       APV->Gain = APV->FitMPV / MPVmean;
0352       if (APV->SubDet > 2)
0353         GOOD++;
0354     } else {
0355       APV->Gain = APV->PreviousGain;
0356       if (APV->SubDet > 2)
0357         BAD++;
0358     }
0359     if (APV->Gain <= 0)
0360       APV->Gain = 1;
0361
0362     delete Proj;
0363   }
0364   printf("\n");
0365 }
0366
0367 //********************************************************************************//
0368 // ------------ method called once each job just before starting event loop  ------------
0369 void SiStripApvGainInspector::checkBookAPVColls(const edm::EventSetup& es) {
0370   tkGeom_ = es.getHandle(tkGeomToken_);
0371   const TrackerGeometry* newBareTkGeomPtr = &(*tkGeom_);
0372   if (newBareTkGeomPtr == bareTkGeomPtr_)
0373     return;  // already filled APVColls, nothing changed
0374
0375   if (!bareTkGeomPtr_) {  // pointer not yet set: called the first time => fill the APVColls
0376     auto const& Det = newBareTkGeomPtr->dets();
0377
0378     unsigned int Index = 0;
0379
0380     for (unsigned int i = 0; i < Det.size(); i++) {
0381       DetId Detid = Det[i]->geographicalId();
0382       int SubDet = Detid.subdetId();
0383
0384       if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
0385           SubDet == StripSubdetector::TEC) {
0386         auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
0387         if (!DetUnit)
0388           continue;
0389
0390         const StripTopology& Topo = DetUnit->specificTopology();
0391         unsigned int NAPV = Topo.nstrips() / 128;
0392
0393         for (unsigned int j = 0; j < NAPV; j++) {
0394           auto APV = std::make_shared<stAPVGain>();
0395           APV->Index = Index;
0396           APV->Bin = -1;
0397           APV->DetId = Detid.rawId();
0398           APV->APVId = j;
0399           APV->SubDet = SubDet;
0400           APV->FitMPV = -1;
0401           APV->FitMPVErr = -1;
0402           APV->FitWidth = -1;
0403           APV->FitWidthErr = -1;
0404           APV->FitChi2 = -1;
0405           APV->FitNorm = -1;
0406           APV->Gain = -1;
0407           APV->PreviousGain = 1;
0408           APV->PreviousGainTick = 1;
0409           APV->x = DetUnit->position().basicVector().x();
0410           APV->y = DetUnit->position().basicVector().y();
0411           APV->z = DetUnit->position().basicVector().z();
0412           APV->Eta = DetUnit->position().basicVector().eta();
0413           APV->Phi = DetUnit->position().basicVector().phi();
0414           APV->R = DetUnit->position().basicVector().transverse();
0415           APV->Thickness = DetUnit->surface().bounds().thickness();
0416           APV->NEntries = 0;
0417           APV->isMasked = false;
0418
0419           APVsCollOrdered.push_back(APV);
0420           APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0421           Index++;
0422           NStripAPVs++;
0423         }  // loop on APVs
0424       }    // if is Strips
0425     }      // loop on dets
0426
0427     for (unsigned int i = 0; i < Det.size();
0428          i++) {  //Make two loop such that the Pixel information is added at the end --> make transition simpler
0429       DetId Detid = Det[i]->geographicalId();
0430       int SubDet = Detid.subdetId();
0431       if (SubDet == PixelSubdetector::PixelBarrel || SubDet == PixelSubdetector::PixelEndcap) {
0432         auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
0433         if (!DetUnit)
0434           continue;
0435
0436         const PixelTopology& Topo = DetUnit->specificTopology();
0437         unsigned int NROCRow = Topo.nrows() / (80.);
0438         unsigned int NROCCol = Topo.ncolumns() / (52.);
0439
0440         for (unsigned int j = 0; j < NROCRow; j++) {
0441           for (unsigned int i = 0; i < NROCCol; i++) {
0442             auto APV = std::make_shared<stAPVGain>();
0443             APV->Index = Index;
0444             APV->Bin = -1;
0445             APV->DetId = Detid.rawId();
0446             APV->APVId = (j << 3 | i);
0447             APV->SubDet = SubDet;
0448             APV->FitMPV = -1;
0449             APV->FitMPVErr = -1;
0450             APV->FitWidth = -1;
0451             APV->FitWidthErr = -1;
0452             APV->FitChi2 = -1;
0453             APV->Gain = -1;
0454             APV->PreviousGain = 1;
0455             APV->PreviousGainTick = 1;
0456             APV->x = DetUnit->position().basicVector().x();
0457             APV->y = DetUnit->position().basicVector().y();
0458             APV->z = DetUnit->position().basicVector().z();
0459             APV->Eta = DetUnit->position().basicVector().eta();
0460             APV->Phi = DetUnit->position().basicVector().phi();
0461             APV->R = DetUnit->position().basicVector().transverse();
0462             APV->Thickness = DetUnit->surface().bounds().thickness();
0463             APV->isMasked = false;  //SiPixelQuality_->IsModuleBad(Detid.rawId());
0464             APV->NEntries = 0;
0465
0466             APVsCollOrdered.push_back(APV);
0467             APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0468             Index++;
0469             NPixelDets++;
0470
0471           }  // loop on ROC cols
0472         }    // loop on ROC rows
0473       }      // if Pixel
0474     }        // loop on Dets
0475   }          //if (!bareTkGeomPtr_) ...
0476   bareTkGeomPtr_ = newBareTkGeomPtr;
0477 }
0478
0479 void SiStripApvGainInspector::storeOnTree(TFileService* tfs) {
0480   unsigned int tree_Index;
0481   unsigned int tree_Bin;
0482   unsigned int tree_DetId;
0483   unsigned char tree_APVId;
0484   unsigned char tree_SubDet;
0485   float tree_x;
0486   float tree_y;
0487   float tree_z;
0488   float tree_Eta;
0489   float tree_R;
0490   float tree_Phi;
0491   float tree_Thickness;
0492   float tree_FitMPV;
0493   float tree_FitMPVErr;
0494   float tree_FitWidth;
0495   float tree_FitWidthErr;
0496   float tree_FitChi2NDF;
0497   float tree_FitNorm;
0498   double tree_Gain;
0499   double tree_PrevGain;
0500   double tree_PrevGainTick;
0501   double tree_NEntries;
0502   bool tree_isMasked;
0503
0504   TTree* MyTree;
0505   MyTree = tfs->make<TTree>("APVGain", "APVGain");
0506   MyTree->Branch("Index", &tree_Index, "Index/i");
0507   MyTree->Branch("Bin", &tree_Bin, "Bin/i");
0508   MyTree->Branch("DetId", &tree_DetId, "DetId/i");
0509   MyTree->Branch("APVId", &tree_APVId, "APVId/b");
0510   MyTree->Branch("SubDet", &tree_SubDet, "SubDet/b");
0511   MyTree->Branch("x", &tree_x, "x/F");
0512   MyTree->Branch("y", &tree_y, "y/F");
0513   MyTree->Branch("z", &tree_z, "z/F");
0514   MyTree->Branch("Eta", &tree_Eta, "Eta/F");
0515   MyTree->Branch("R", &tree_R, "R/F");
0516   MyTree->Branch("Phi", &tree_Phi, "Phi/F");
0517   MyTree->Branch("Thickness", &tree_Thickness, "Thickness/F");
0518   MyTree->Branch("FitMPV", &tree_FitMPV, "FitMPV/F");
0519   MyTree->Branch("FitMPVErr", &tree_FitMPVErr, "FitMPVErr/F");
0520   MyTree->Branch("FitWidth", &tree_FitWidth, "FitWidth/F");
0521   MyTree->Branch("FitWidthErr", &tree_FitWidthErr, "FitWidthErr/F");
0522   MyTree->Branch("FitChi2NDF", &tree_FitChi2NDF, "FitChi2NDF/F");
0523   MyTree->Branch("FitNorm", &tree_FitNorm, "FitNorm/F");
0524   MyTree->Branch("Gain", &tree_Gain, "Gain/D");
0525   MyTree->Branch("PrevGain", &tree_PrevGain, "PrevGain/D");
0526   MyTree->Branch("PrevGainTick", &tree_PrevGainTick, "PrevGainTick/D");
0527   MyTree->Branch("NEntries", &tree_NEntries, "NEntries/D");
0528   MyTree->Branch("isMasked", &tree_isMasked, "isMasked/O");
0529
0530   uint32_t cachedId(0);
0531   SiStripMiscalibrate::Entry gain_ratio;
0532   SiStripMiscalibrate::Entry o_gain;
0533   SiStripMiscalibrate::Entry n_gain;
0534   SiStripMiscalibrate::Entry mpv;
0535   SiStripMiscalibrate::Entry mpv_err;
0536   SiStripMiscalibrate::Entry entries;
0537   SiStripMiscalibrate::Entry fitChi2;
0538
0539   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0540     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0541     if (APV == nullptr)
0542       continue;
0543     //     printf(      "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
0544     //fprintf(Gains,"%i | %i | PreviousGain = %7.5f(tick) x %7.5f(particle) NewGain (particle) = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain,APV->Gain, APV->NEntries);
0545
0546     // do not fill the Pixel
0547     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0548       continue;
0549
0550     tree_Index = APV->Index;
0551     tree_Bin = Charge_Vs_Index->GetXaxis()->FindBin(APV->Index);
0552     tree_DetId = APV->DetId;
0553     tree_APVId = APV->APVId;
0554     tree_SubDet = APV->SubDet;
0555     tree_x = APV->x;
0556     tree_y = APV->y;
0557     tree_z = APV->z;
0558     tree_Eta = APV->Eta;
0559     tree_R = APV->R;
0560     tree_Phi = APV->Phi;
0561     tree_Thickness = APV->Thickness;
0562     tree_FitMPV = APV->FitMPV;
0563     tree_FitMPVErr = APV->FitMPVErr;
0564     tree_FitWidth = APV->FitWidth;
0565     tree_FitWidthErr = APV->FitWidthErr;
0566     tree_FitChi2NDF = APV->FitChi2;
0567     tree_FitNorm = APV->FitNorm;
0568     tree_Gain = APV->Gain;
0569     tree_PrevGain = APV->PreviousGain;
0570     tree_PrevGainTick = APV->PreviousGainTick;
0571     tree_NEntries = APV->NEntries;
0572     tree_isMasked = APV->isMasked;
0573
0574     // flush the counters
0575     if (cachedId != 0 && tree_DetId != cachedId) {
0576       //ratio_map->fill(cachedId,gain_ratio.mean());
0577       ratio_map->fill(cachedId, o_gain.mean() / n_gain.mean());
0578       old_payload_map->fill(cachedId, o_gain.mean());
0579       new_payload_map->fill(cachedId, n_gain.mean());
0580
0581       if (entries.mean() > 0) {
0582         mpv_map->fill(cachedId, mpv.mean());
0583         mpv_err_map->fill(cachedId, mpv_err.mean());
0584         entries_map->fill(cachedId, log10(entries.mean()));
0585         if (fitChi2.mean() > 0) {
0586           fitChi2_map->fill(cachedId, log10(fitChi2.mean()));
0587         } else {
0588           fitChi2_map->fill(cachedId, -1);
0589         }
0590       }
0591
0592       gain_ratio.reset();
0593       o_gain.reset();
0594       n_gain.reset();
0595
0596       mpv.reset();
0597       mpv_err.reset();
0598       entries.reset();
0599       fitChi2.reset();
0600     }
0601
0602     cachedId = tree_DetId;
0603     gain_ratio.add(tree_PrevGain / tree_Gain);
0604     o_gain.add(tree_PrevGain);
0605     n_gain.add(tree_Gain);
0606     mpv.add(tree_FitMPV);
0607     mpv_err.add(tree_FitMPVErr);
0608     entries.add(tree_NEntries);
0609     fitChi2.add(tree_FitChi2NDF);
0610
0611     if (tree_DetId == 402673324) {
0612       printf("%i | %i : %f --> %f  (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
0613     }
0614
0615     MyTree->Fill();
0616   }
0617 }
0618
0619 //********************************************************************************//
0620 void SiStripApvGainInspector::checkAndRetrieveTopology(const edm::EventSetup& setup) {
0621   if (!tTopo_) {
0622     edm::ESHandle<TrackerTopology> TopoHandle = setup.getHandle(tTopoToken_);
0623     tTopo_ = TopoHandle.product();
0624   }
0625 }
0626
0627 //********************************************************************************//
0628 void SiStripApvGainInspector::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange) {
0629   FitResults[0] = -0.5;  //MPV
0630   FitResults[1] = 0;     //MPV error
0631   FitResults[2] = -0.5;  //Width
0632   FitResults[3] = 0;     //Width error
0633   FitResults[4] = -0.5;  //Fit Chi2/NDF
0634   FitResults[5] = 0;     //Normalization
0635
0636   if (InputHisto->GetEntries() < minNrEntries)
0637     return;
0638
0639   // perform fit with standard landau
0640   TF1 MyLandau("MyLandau", "landau", LowRange, HighRange);
0641   MyLandau.SetParameter(1, 300);
0642   InputHisto->Fit(&MyLandau, "QR WW");
0643
0644   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0645   FitResults[0] = MyLandau.GetParameter(1);                     //MPV
0646   FitResults[1] = MyLandau.GetParError(1);                      //MPV error
0647   FitResults[2] = MyLandau.GetParameter(2);                     //Width
0648   FitResults[3] = MyLandau.GetParError(2);                      //Width error
0649   FitResults[4] = MyLandau.GetChisquare() / MyLandau.GetNDF();  //Fit Chi2/NDF
0650   FitResults[5] = MyLandau.GetParameter(0);
0651 }
0652
0653 void SiStripApvGainInspector::doFakeFit(TH1* InputHisto, double* FitResults) {
0654   FitResults[0] = -0.5;  //MPV
0655   FitResults[1] = 0;     //MPV error
0656   FitResults[2] = -0.5;  //Width
0657   FitResults[3] = 0;     //Width error
0658   FitResults[4] = -0.5;  //Fit Chi2/NDF
0659   FitResults[5] = 0;     //Normalization
0660 }
0661
0662 //********************************************************************************//
0663 double SiStripApvGainInspector::langaufun(Double_t* x, Double_t* par)
0664 //********************************************************************************//
0665 {
0666   //Fit parameters:
0667   //par[0]=Width (scale) parameter of Landau density
0668   //par[1]=Most Probable (MP, location) parameter of Landau density
0669   //par[2]=Total area (integral -inf to inf, normalization constant)
0670   //par[3]=Width (sigma) of convoluted Gaussian function
0671   //
0672   //In the Landau distribution (represented by the CERNLIB approximation),
0673   //the maximum is located at x=-0.22278298 with the location parameter=0.
0674   //This shift is corrected within this function, so that the actual
0675   //maximum is identical to the MP parameter.
0676
0677   // Numeric constants
0678   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0679   Double_t mpshift = -0.22278298;       // Landau maximum location
0680
0681   // Control constants
0682   Double_t np = 100.0;  // number of convolution steps
0683   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0684
0685   // Variables
0686   Double_t xx;
0687   Double_t mpc;
0688   Double_t fland;
0689   Double_t sum = 0.0;
0690   Double_t xlow, xupp;
0691   Double_t step;
0692   Double_t i;
0693
0694   // MP shift correction
0695   mpc = par[1] - mpshift * par[0];
0696
0697   // Range of convolution integral
0698   xlow = x[0] - sc * par[3];
0699   xupp = x[0] + sc * par[3];
0700
0701   step = (xupp - xlow) / np;
0702
0703   // Convolution integral of Landau and Gaussian by sum
0704   for (i = 1.0; i <= np / 2; i++) {
0705     xx = xlow + (i - .5) * step;
0706     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0707     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0708
0709     xx = xupp - (i - .5) * step;
0710     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0711     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0712   }
0713
0714   return (par[2] * step * sum * invsq2pi / par[3]);
0715 }
0716
0717 //********************************************************************************//
0718 void SiStripApvGainInspector::getPeakOfLanGau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange) {
0719   FitResults[0] = -0.5;  //MPV
0720   FitResults[1] = 0;     //MPV error
0721   FitResults[2] = -0.5;  //Width
0722   FitResults[3] = 0;     //Width error
0723   FitResults[4] = -0.5;  //Fit Chi2/NDF
0724   FitResults[5] = 0;     //Normalization
0725
0726   if (InputHisto->GetEntries() < minNrEntries)
0727     return;
0728
0729   // perform fit with standard landau
0730   TF1 MyLandau("MyLandau", "landau", LowRange, HighRange);
0731   MyLandau.SetParameter(1, 300);
0732   InputHisto->Fit(&MyLandau, "QR WW");
0733
0734   double startvalues[4] = {100, 300, 10000, 100};
0735   double parlimitslo[4] = {0, 250, 10, 0};
0736   double parlimitshi[4] = {200, 350, 1000000, 200};
0737
0738   TF1 MyLangau("MyLanGau", langaufun, LowRange, HighRange, 4);
0739
0740   MyLangau.SetParameters(startvalues);
0741   MyLangau.SetParNames("Width", "MP", "Area", "GSigma");
0742
0743   for (unsigned int i = 0; i < 4; i++) {
0744     MyLangau.SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0745   }
0746
0747   InputHisto->Fit("MyLanGau", "QRB0");  // fit within specified range, use ParLimits, do not plot
0748
0749   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0750   FitResults[0] = MyLangau.GetParameter(1);                     //MPV
0751   FitResults[1] = MyLangau.GetParError(1);                      //MPV error
0752   FitResults[2] = MyLangau.GetParameter(0);                     //Width
0753   FitResults[3] = MyLangau.GetParError(0);                      //Width error
0754   FitResults[4] = MyLangau.GetChisquare() / MyLangau.GetNDF();  //Fit Chi2/NDF
0755   FitResults[5] = MyLangau.GetParameter(3);
0756 }
0757
0758 //********************************************************************************//
0759 void SiStripApvGainInspector::getPeakOfLandauAroundMax(TH1* InputHisto,
0760                                                        double* FitResults,
0761                                                        double LowRange,
0762                                                        double HighRange) {
0763   FitResults[0] = -0.5;  //MPV
0764   FitResults[1] = 0;     //MPV error
0765   FitResults[2] = -0.5;  //Width
0766   FitResults[3] = 0;     //Width error
0767   FitResults[4] = -0.5;  //Fit Chi2/NDF
0768   FitResults[5] = 0;     //Normalization
0769
0770   if (InputHisto->GetEntries() < minNrEntries)
0771     return;
0772
0773   int maxbin = InputHisto->GetMaximumBin();
0774   int maxbin2 = -9999.;
0775
0776   if (InputHisto->GetBinContent(maxbin - 1) > InputHisto->GetBinContent(maxbin + 1)) {
0777     maxbin2 = maxbin - 1;
0778   } else {
0779     maxbin2 = maxbin + 1;
0780   }
0781
0782   float maxbincenter = (InputHisto->GetBinCenter(maxbin) + InputHisto->GetBinCenter(maxbin2)) / 2;
0783
0784   TF1 MyLandau("MyLandau", "[2]*TMath::Landau(x,[0],[1],0)", maxbincenter - LowRange, maxbincenter + HighRange);
0785   // TF1 MyLandau("MyLandau", "landau", LowRange, HighRange);
0786   // MyLandau.SetParameter(1, 300);
0787   InputHisto->Fit(&MyLandau, "QR WW");
0788
0789   MyLandau.SetParameter(0, maxbincenter);
0790   MyLandau.SetParameter(1, maxbincenter / 10.);
0791   MyLandau.SetParameter(2, InputHisto->GetMaximum());
0792
0793   float mpv = MyLandau.GetParameter(1);
0794   MyLandau.SetParameter(1, mpv);
0795   //InputHisto->Rebin(3);
0796   InputHisto->Fit(&MyLandau, "QOR", "", mpv - 50, mpv + 100);
0797
0798   InputHisto->Fit(&MyLandau, "QOR", "", maxbincenter - LowRange, maxbincenter + HighRange);
0799   InputHisto->Fit(&MyLandau, "QOR", "", maxbincenter - LowRange, maxbincenter + HighRange);
0800
0801   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0802   FitResults[0] = MyLandau.GetParameter(1);                     //MPV
0803   FitResults[1] = MyLandau.GetParError(1);                      //MPV error
0804   FitResults[2] = MyLandau.GetParameter(2);                     //Width
0805   FitResults[3] = MyLandau.GetParError(2);                      //Width error
0806   FitResults[4] = MyLandau.GetChisquare() / MyLandau.GetNDF();  //Fit Chi2/NDF
0807   FitResults[5] = MyLandau.GetParameter(0);
0808 }
0809
0810 //********************************************************************************//
0811 bool SiStripApvGainInspector::isGoodLandauFit(double* FitResults) {
0812   if (FitResults[0] <= 0)
0813     return false;
0814   //   if(FitResults[1] > MaxMPVError   )return false;
0815   //   if(FitResults[4] > MaxChi2OverNDF)return false;
0816   return true;
0817 }
0818
0819 /*--------------------------------------------------------------------*/
0820 void SiStripApvGainInspector::makeNicePlotStyle(TH1F* plot)
0821 /*--------------------------------------------------------------------*/
0822 {
0823   plot->GetXaxis()->CenterTitle(true);
0824   plot->GetYaxis()->CenterTitle(true);
0825   plot->GetXaxis()->SetTitleFont(42);
0826   plot->GetYaxis()->SetTitleFont(42);
0827   plot->GetXaxis()->SetTitleSize(0.05);
0828   plot->GetYaxis()->SetTitleSize(0.05);
0829   plot->GetXaxis()->SetTitleOffset(0.9);
0830   plot->GetYaxis()->SetTitleOffset(1.3);
0831   plot->GetXaxis()->SetLabelFont(42);
0832   plot->GetYaxis()->SetLabelFont(42);
0833   plot->GetYaxis()->SetLabelSize(.05);
0834   plot->GetXaxis()->SetLabelSize(.05);
0835 }
0836
0837 //********************************************************************************//
0838 std::unique_ptr<SiStripApvGain> SiStripApvGainInspector::getNewObject() {
0839   std::unique_ptr<SiStripApvGain> obj = std::make_unique<SiStripApvGain>();
0840
0841   std::vector<float> theSiStripVector;
0842   unsigned int PreviousDetId = 0;
0843   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0844     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0845     if (APV == nullptr) {
0846       printf("Bug\n");
0847       continue;
0848     }
0849     if (APV->SubDet <= 2)
0850       continue;
0851     if (APV->DetId != PreviousDetId) {
0852       if (!theSiStripVector.empty()) {
0853         SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0854         if (!obj->put(PreviousDetId, range))
0855           printf("Bug to put detId = %i\n", PreviousDetId);
0856       }
0857       theSiStripVector.clear();
0858       PreviousDetId = APV->DetId;
0859     }
0860     theSiStripVector.push_back(APV->Gain);
0861
0862     LogDebug("SiStripGainsPCLHarvester") << " DetId: " << APV->DetId << " APV:   " << APV->APVId
0863                                          << " Gain:  " << APV->Gain << std::endl;
0864   }
0865   if (!theSiStripVector.empty()) {
0866     SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0867     if (!obj->put(PreviousDetId, range))
0868       printf("Bug to put detId = %i\n", PreviousDetId);
0869   }
0870
0871   return obj;
0872 }
0873
0874 // ------------ method called once each job just before starting event loop  ------------
0875 void SiStripApvGainInspector::beginJob() {
0876   TFileDirectory control_dir = tfs->mkdir("Control");
0877   //DA.cd();
0878   hControl = this->bookQualityMonitor(control_dir);
0879 }
0880
0881 // ------------ method called once each job just after ending the event loop  ------------
0882 void SiStripApvGainInspector::endJob() {
0883   edm::LogVerbatim("SelectedModules") << "Selected APVs:" << histoMap_.size() << std::endl;
0884   for (const auto& plot : histoMap_) {
0885     TCanvas* c1 = new TCanvas(Form("c1_%i_%i", plot.first.second, plot.first.first),
0886                               Form("c1_%i_%i", plot.first.second, plot.first.first),
0887                               800,
0888                               600);
0889     // Define common things for the different fits
0890
0891     gStyle->SetOptFit(1011);
0892     c1->Clear();
0893
0894     c1->SetLeftMargin(0.15);
0895     c1->SetRightMargin(0.10);
0896     plot.second->SetTitle(Form("Cluster Charge (%i,%i)", plot.first.second, plot.first.first));
0897     plot.second->GetXaxis()->SetTitle("Normalized Cluster Charge [ADC counts/mm]");
0898     plot.second->GetYaxis()->SetTitle("On-track clusters");
0899     plot.second->GetXaxis()->SetRangeUser(0., 1000.);
0900
0901     this->makeNicePlotStyle(plot.second);
0902     plot.second->Draw();
0903     edm::LogVerbatim("SelectedModules") << " DetId: " << plot.first.second << " (" << plot.first.first << ")"
0904                                         << std::endl;
0905     ;
0906
0907     c1->Print(Form("c1_%i_%i.png", plot.first.second, plot.first.first));
0908     c1->Print(Form("c1_%i_%i.pdf", plot.first.second, plot.first.first));
0909   }
0910
0911   tfs = edm::Service<TFileService>().operator->();
0912   storeOnTree(tfs);
0913
0914   auto range = SiStripMiscalibrate::getTruncatedRange(ratio_map.get());
0915
0916   ratio_map->save(true, range.first, range.second, "G2_gain_ratio_map.pdf");
0917   ratio_map->save(true, range.first, range.second, "G2_gain_ratio_map.png");
0918
0919   range = SiStripMiscalibrate::getTruncatedRange(old_payload_map.get());
0920
0921   old_payload_map->save(true, range.first, range.second, "starting_G2_gain_payload_map.pdf");
0922   old_payload_map->save(true, range.first, range.second, "starting_G2_gain_payload_map.png");
0923
0924   range = SiStripMiscalibrate::getTruncatedRange(new_payload_map.get());
0925
0926   new_payload_map->save(true, range.first, range.second, "new_G2_gain_payload_map.pdf");
0927   new_payload_map->save(true, range.first, range.second, "new_G2_gain_payload_map.png");
0928
0929   mpv_map->save(true, 250, 350., "mpv_map.pdf");
0930   mpv_map->save(true, 250, 350., "mpv_map.png");
0931
0932   mpv_err_map->save(true, 0., 3., "mpv_err_map.pdf");
0933   mpv_err_map->save(true, 0., 3., "mpv_err_map.png");
0934
0935   entries_map->save(true, 0, 0, "entries_map.pdf");
0936   entries_map->save(true, 0, 0, "entries_map.png");
0937
0938   fitChi2_map->save(true, 0., 0., "fitChi2_map.pdf");
0939   fitChi2_map->save(true, 0., 0., "fitChi2_map.png");
0940
0941   fillQualityMonitor();
0942
0943   std::unique_ptr<SiStripApvGain> theAPVGains = this->getNewObject();
0944
0945   // write out the APVGains record
0946   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0947
0948   if (poolDbService.isAvailable())
0949     poolDbService->writeOneIOV(theAPVGains.get(), poolDbService->currentTime(), "SiStripApvGainRcd");
0950   else
0951     throw std::runtime_error("PoolDBService required.");
0952 }
0953
0954 std::map<std::string, TH1*> SiStripApvGainInspector::bookQualityMonitor(const TFileDirectory& dir) {
0955   int MPVbin = 300;
0956   float MPVmin = 0.;
0957   float MPVmax = 600.;
0958
0959   TH1F::SetDefaultSumw2(kTRUE);
0960   std::map<std::string, TH1*> h;
0961
0962   h["MPV_Vs_EtaTIB"] = dir.make<TH2F>("MPVvsEtaTIB", "MPV vs Eta TIB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0963   h["MPV_Vs_EtaTID"] = dir.make<TH2F>("MPVvsEtaTID", "MPV vs Eta TID", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0964   h["MPV_Vs_EtaTOB"] = dir.make<TH2F>("MPVvsEtaTOB", "MPV vs Eta TOB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0965   h["MPV_Vs_EtaTEC"] = dir.make<TH2F>("MPVvsEtaTEC", "MPV vs Eta TEC", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0966   h["MPV_Vs_EtaTECthin"] = dir.make<TH2F>("MPVvsEtaTEC1", "MPV vs Eta TEC-thin", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0967   h["MPV_Vs_EtaTECthick"] =
0968       dir.make<TH2F>("MPVvsEtaTEC2", "MPV vs Eta TEC-thick", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0969
0970   h["MPV_Vs_PhiTIB"] = dir.make<TH2F>("MPVvsPhiTIB", "MPV vs Phi TIB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0971   h["MPV_Vs_PhiTID"] = dir.make<TH2F>("MPVvsPhiTID", "MPV vs Phi TID", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0972   h["MPV_Vs_PhiTOB"] = dir.make<TH2F>("MPVvsPhiTOB", "MPV vs Phi TOB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0973   h["MPV_Vs_PhiTEC"] = dir.make<TH2F>("MPVvsPhiTEC", "MPV vs Phi TEC", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0974   h["MPV_Vs_PhiTECthin"] =
0975       dir.make<TH2F>("MPVvsPhiTEC1", "MPV vs Phi TEC-thin ", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0976   h["MPV_Vs_PhiTECthick"] =
0977       dir.make<TH2F>("MPVvsPhiTEC2", "MPV vs Phi TEC-thick", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0978
0979   h["NoMPVfit"] = dir.make<TH2F>("NoMPVfit", "Modules with bad Landau Fit", 350, -350, 350, 240, 0, 120);
0980   h["NoMPVmasked"] = dir.make<TH2F>("NoMPVmasked", "Masked Modules", 350, -350, 350, 240, 0, 120);
0981
0982   h["Gains"] = dir.make<TH1F>("Gains", "Gains", 300, 0, 2);
0983   h["MPVs"] = dir.make<TH1F>("MPVs", "MPVs", MPVbin, MPVmin, MPVmax);
0984   h["MPVs320"] = dir.make<TH1F>("MPV_320", "MPV 320 thickness", MPVbin, MPVmin, MPVmax);
0985   h["MPVs500"] = dir.make<TH1F>("MPV_500", "MPV 500 thickness", MPVbin, MPVmin, MPVmax);
0986   h["MPVsTIB"] = dir.make<TH1F>("MPV_TIB", "MPV TIB", MPVbin, MPVmin, MPVmax);
0987   h["MPVsTID"] = dir.make<TH1F>("MPV_TID", "MPV TID", MPVbin, MPVmin, MPVmax);
0988   h["MPVsTIDP"] = dir.make<TH1F>("MPV_TIDP", "MPV TIDP", MPVbin, MPVmin, MPVmax);
0989   h["MPVsTIDM"] = dir.make<TH1F>("MPV_TIDM", "MPV TIDM", MPVbin, MPVmin, MPVmax);
0990   h["MPVsTOB"] = dir.make<TH1F>("MPV_TOB", "MPV TOB", MPVbin, MPVmin, MPVmax);
0991   h["MPVsTEC"] = dir.make<TH1F>("MPV_TEC", "MPV TEC", MPVbin, MPVmin, MPVmax);
0992   h["MPVsTECP"] = dir.make<TH1F>("MPV_TECP", "MPV TECP", MPVbin, MPVmin, MPVmax);
0993   h["MPVsTECM"] = dir.make<TH1F>("MPV_TECM", "MPV TECM", MPVbin, MPVmin, MPVmax);
0994   h["MPVsTECthin"] = dir.make<TH1F>("MPV_TEC1", "MPV TEC thin", MPVbin, MPVmin, MPVmax);
0995   h["MPVsTECthick"] = dir.make<TH1F>("MPV_TEC2", "MPV TEC thick", MPVbin, MPVmin, MPVmax);
0996   h["MPVsTECP1"] = dir.make<TH1F>("MPV_TECP1", "MPV TECP thin ", MPVbin, MPVmin, MPVmax);
0997   h["MPVsTECP2"] = dir.make<TH1F>("MPV_TECP2", "MPV TECP thick", MPVbin, MPVmin, MPVmax);
0998   h["MPVsTECM1"] = dir.make<TH1F>("MPV_TECM1", "MPV TECM thin", MPVbin, MPVmin, MPVmax);
0999   h["MPVsTECM2"] = dir.make<TH1F>("MPV_TECM2", "MPV TECM thick", MPVbin, MPVmin, MPVmax);
1000
1001   h["MPVError"] = dir.make<TH1F>("MPVError", "MPV Error", 150, 0, 150);
1002   h["MPVErrorVsMPV"] = dir.make<TH2F>("MPVErrorVsMPV", "MPV Error vs MPV", 300, 0, 600, 150, 0, 150);
1003   h["MPVErrorVsEta"] = dir.make<TH2F>("MPVErrorVsEta", "MPV Error vs Eta", 50, -3.0, 3.0, 150, 0, 150);
1004   h["MPVErrorVsPhi"] = dir.make<TH2F>("MPVErrorVsPhi", "MPV Error vs Phi", 50, -3.4, 3.4, 150, 0, 150);
1005   h["MPVErrorVsN"] = dir.make<TH2F>("MPVErrorVsN", "MPV Error vs N", 500, 0, 1000, 150, 0, 150);
1006
1007   h["DiffWRTPrevGainTIB"] = dir.make<TH1F>("DiffWRTPrevGainTIB", "Diff w.r.t. PrevGain TIB", 250, 0.5, 1.5);
1008   h["DiffWRTPrevGainTID"] = dir.make<TH1F>("DiffWRTPrevGainTID", "Diff w.r.t. PrevGain TID", 250, 0.5, 1.5);
1009   h["DiffWRTPrevGainTOB"] = dir.make<TH1F>("DiffWRTPrevGainTOB", "Diff w.r.t. PrevGain TOB", 250, 0.5, 1.5);
1010   h["DiffWRTPrevGainTEC"] = dir.make<TH1F>("DiffWRTPrevGainTEC", "Diff w.r.t. PrevGain TEC", 250, 0.5, 1.5);
1011
1012   h["GainVsPrevGainTIB"] = dir.make<TH2F>("GainVsPrevGainTIB", "Gain vs PrevGain TIB", 100, 0, 2, 100, 0, 2);
1013   h["GainVsPrevGainTID"] = dir.make<TH2F>("GainVsPrevGainTID", "Gain vs PrevGain TID", 100, 0, 2, 100, 0, 2);
1014   h["GainVsPrevGainTOB"] = dir.make<TH2F>("GainVsPrevGainTOB", "Gain vs PrevGain TOB", 100, 0, 2, 100, 0, 2);
1015   h["GainVsPrevGainTEC"] = dir.make<TH2F>("GainVsPrevGainTEC", "Gain vs PrevGain TEC", 100, 0, 2, 100, 0, 2);
1016
1017   return h;
1018 }
1019
1020 void SiStripApvGainInspector::fillQualityMonitor() {
1021   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
1022     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
1023     if (APV == nullptr)
1024       continue;
1025
1026     //unsigned int Index = APV->Index;
1027     //unsigned int DetId = APV->DetId;
1028     unsigned int SubDet = APV->SubDet;
1029     float z = APV->z;
1030     float Eta = APV->Eta;
1031     float R = APV->R;
1032     float Phi = APV->Phi;
1033     float Thickness = APV->Thickness;
1034     double FitMPV = APV->FitMPV;
1035     double FitMPVErr = APV->FitMPVErr;
1036     double Gain = APV->Gain;
1037     double NEntries = APV->NEntries;
1038     double PreviousGain = APV->PreviousGain;
1039
1040     if (SubDet < 3)
1041       continue;  // avoid to loop over Pixel det id
1042
1043     if (FitMPV <= 0.) {  // No fit of MPV
1044       if (APV->isMasked)
1045         fill2D(hControl, "NoMPVmasked", z, R);
1046       else
1047         fill2D(hControl, "NoMPVfit", z, R);
1048     } else {  // Fit of MPV
1049       if (FitMPV > 0.)
1050         fill1D(hControl, "Gains", Gain);
1051
1052       fill1D(hControl, "MPVs", FitMPV);
1053       if (Thickness < 0.04)
1054         fill1D(hControl, "MPVs320", FitMPV);
1055       if (Thickness > 0.04)
1056         fill1D(hControl, "MPVs500", FitMPV);
1057
1058       fill1D(hControl, "MPVError", FitMPVErr);
1059       fill2D(hControl, "MPVErrorVsMPV", FitMPV, FitMPVErr);
1060       fill2D(hControl, "MPVErrorVsEta", Eta, FitMPVErr);
1061       fill2D(hControl, "MPVErrorVsPhi", Phi, FitMPVErr);
1062       fill2D(hControl, "MPVErrorVsN", NEntries, FitMPVErr);
1063
1064       if (SubDet == 3) {
1065         fill2D(hControl, "MPV_Vs_EtaTIB", Eta, FitMPV);
1066         fill2D(hControl, "MPV_Vs_PhiTIB", Phi, FitMPV);
1067         fill1D(hControl, "MPVsTIB", FitMPV);
1068
1069       } else if (SubDet == 4) {
1070         fill2D(hControl, "MPV_Vs_EtaTID", Eta, FitMPV);
1071         fill2D(hControl, "MPV_Vs_PhiTID", Phi, FitMPV);
1072         fill1D(hControl, "MPVsTID", FitMPV);
1073         if (Eta < 0.)
1074           fill1D(hControl, "MPVsTIDM", FitMPV);
1075         if (Eta > 0.)
1076           fill1D(hControl, "MPVsTIDP", FitMPV);
1077
1078       } else if (SubDet == 5) {
1079         fill2D(hControl, "MPV_Vs_EtaTOB", Eta, FitMPV);
1080         fill2D(hControl, "MPV_Vs_PhiTOB", Phi, FitMPV);
1081         fill1D(hControl, "MPVsTOB", FitMPV);
1082
1083       } else if (SubDet == 6) {
1084         fill2D(hControl, "MPV_Vs_EtaTEC", Eta, FitMPV);
1085         fill2D(hControl, "MPV_Vs_PhiTEC", Phi, FitMPV);
1086         fill1D(hControl, "MPVsTEC", FitMPV);
1087         if (Eta < 0.)
1088           fill1D(hControl, "MPVsTECM", FitMPV);
1089         if (Eta > 0.)
1090           fill1D(hControl, "MPVsTECP", FitMPV);
1091         if (Thickness < 0.04) {
1092           fill2D(hControl, "MPV_Vs_EtaTECthin", Eta, FitMPV);
1093           fill2D(hControl, "MPV_Vs_PhiTECthin", Phi, FitMPV);
1094           fill1D(hControl, "MPVsTECthin", FitMPV);
1095           if (Eta > 0.)
1096             fill1D(hControl, "MPVsTECP1", FitMPV);
1097           if (Eta < 0.)
1098             fill1D(hControl, "MPVsTECM1", FitMPV);
1099         }
1100         if (Thickness > 0.04) {
1101           fill2D(hControl, "MPV_Vs_EtaTECthick", Eta, FitMPV);
1102           fill2D(hControl, "MPV_Vs_PhiTECthick", Phi, FitMPV);
1103           fill1D(hControl, "MPVsTECthick", FitMPV);
1104           if (Eta > 0.)
1105             fill1D(hControl, "MPVsTECP2", FitMPV);
1106           if (Eta < 0.)
1107             fill1D(hControl, "MPVsTECM2", FitMPV);
1108         }
1109       }
1110     }
1111
1112     if (SubDet == 3 && PreviousGain != 0.)
1113       fill1D(hControl, "DiffWRTPrevGainTIB", Gain / PreviousGain);
1114     else if (SubDet == 4 && PreviousGain != 0.)
1115       fill1D(hControl, "DiffWRTPrevGainTID", Gain / PreviousGain);
1116     else if (SubDet == 5 && PreviousGain != 0.)
1117       fill1D(hControl, "DiffWRTPrevGainTOB", Gain / PreviousGain);
1118     else if (SubDet == 6 && PreviousGain != 0.)
1119       fill1D(hControl, "DiffWRTPrevGainTEC", Gain / PreviousGain);
1120
1121     if (SubDet == 3)
1122       fill2D(hControl, "GainVsPrevGainTIB", PreviousGain, Gain);
1123     else if (SubDet == 4)
1124       fill2D(hControl, "GainVsPrevGainTID", PreviousGain, Gain);
1125     else if (SubDet == 5)
1126       fill2D(hControl, "GainVsPrevGainTOB", PreviousGain, Gain);
1127     else if (SubDet == 6)
1128       fill2D(hControl, "GainVsPrevGainTEC", PreviousGain, Gain);
1129
1130   }  // loop on the APV collections
1131 }
1132
1133 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
1134 void SiStripApvGainInspector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1135   edm::ParameterSetDescription desc;
1136   desc.addUntracked<std::string>("inputFile", {});
1137   desc.addUntracked<double>("minNrEntries", 20);
1138   desc.add<int>("fitMode", 2)
1139       ->setComment("fit mode. Available options: 1: landau\n 2: landau around max\n 3:landau&gaus convo\n 4: fake");
1140   desc.addUntracked<std::vector<unsigned int>>("selectedModules", {});
1141   descriptions.addWithDefaultLabel(desc);
1142 }
1143
1144 //define this as a plug-in
1145 DEFINE_FWK_MODULE(SiStripApvGainInspector);