Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-24 01:11:13

0001 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0002 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0003 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
0004 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0005 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
0006 #include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h"
0007 #include "DQMServices/Core/interface/DQMStore.h"
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0010 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/EDGetToken.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0022 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0023 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0024 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0025 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0026 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 
0029 #include <TChain.h>
0030 
0031 /// user includes
0032 #include "CalibTracker/SiStripChannelGain/interface/APVGainStruct.h"
0033 #include "CalibTracker/SiStripChannelGain/interface/APVGainHelpers.h"
0034 
0035 //
0036 /**\class SiStripGainsCalibTreeWorker SiStripGainsCalibTreeWorker.cc
0037    Description: Fill DQM histograms with SiStrip Charge normalized to path length
0038 */
0039 //
0040 //  Original Author: L. Quertermont (calibration algorithm)
0041 //  Contributors:    M. Verzetti    (data access)
0042 //                   A. Di Mattia   (PCL multi stream processing and monitoring)
0043 //                   M. Delcourt    (monitoring)
0044 //                   M. Musich      (migration to thread-safe DQMStore access)
0045 //                   P. David       (adapted to calibtrees)
0046 //
0047 //  Created:  Wed, 12 Apr 2017 14:46:48 GMT
0048 //
0049 
0050 class SiStripGainsCalibTreeWorker : public DQMGlobalEDAnalyzer<APVGain::APVGainHistograms> {
0051 public:
0052   explicit SiStripGainsCalibTreeWorker(const edm::ParameterSet&);
0053 
0054   void bookHistograms(DQMStore::IBooker&,
0055                       edm::Run const&,
0056                       edm::EventSetup const&,
0057                       APVGain::APVGainHistograms&) const override;
0058   void dqmAnalyze(edm::Event const&, edm::EventSetup const&, APVGain::APVGainHistograms const&) const override;
0059 
0060   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0061 
0062 private:
0063   void dqmBeginRun(edm::Run const&, edm::EventSetup const&, APVGain::APVGainHistograms&) const override;
0064   void checkBookAPVColls(const TrackerGeometry* bareTkGeomPtr, APVGain::APVGainHistograms& histograms) const;
0065 
0066   std::vector<std::string> dqm_tag_;
0067 
0068   int statCollectionFromMode(const char* tag) const;
0069 
0070   double MinTrackMomentum;
0071   double MaxTrackMomentum;
0072   double MinTrackEta;
0073   double MaxTrackEta;
0074   unsigned int MaxNrStrips;
0075   unsigned int MinTrackHits;
0076   double MaxTrackChiOverNdf;
0077   int MaxTrackingIteration;
0078   bool AllowSaturation;
0079   bool FirstSetOfConstants;
0080   bool Validation;
0081   bool OldGainRemoving;
0082   bool useCalibration;
0083   bool doChargeMonitorPerPlane; /*!< Charge monitor per detector plane */
0084 
0085   mutable bool hasProcessed_;
0086 
0087   std::string m_DQMdir;                  /*!< DQM folder hosting the charge statistics and the monitor plots */
0088   std::string m_calibrationMode;         /*!< Type of statistics for the calibration */
0089   std::vector<std::string> VChargeHisto; /*!< Charge monitor plots to be output */
0090 
0091   // maps histograms index to topology
0092   std::map<unsigned int, APVloc> theTopologyMap;
0093 
0094   std::string EventPrefix_;  //("");
0095   std::string EventSuffix_;  //("");
0096   std::string TrackPrefix_;  //("track");
0097   std::string TrackSuffix_;  //("");
0098   std::string CalibPrefix_;  //("GainCalibration");
0099   std::string CalibSuffix_;  //("");
0100   std::string calibTreeName_;
0101   std::vector<std::string> calibTreeFileNames_;
0102 
0103   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0104   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomTokenBR_, tkGeomToken_;
0105   edm::ESGetToken<SiStripGain, SiStripGainRcd> gainToken_;
0106   edm::ESGetToken<SiStripQuality, SiStripQualityRcd> qualityToken_;
0107 };
0108 
0109 inline int SiStripGainsCalibTreeWorker::statCollectionFromMode(const char* tag) const {
0110   std::vector<std::string>::const_iterator it = dqm_tag_.begin();
0111   while (it != dqm_tag_.end()) {
0112     if (*it == std::string(tag))
0113       return it - dqm_tag_.begin();
0114     it++;
0115   }
0116 
0117   if (std::string(tag).empty())
0118     return 0;  // return StdBunch calibration mode for backward compatibility
0119 
0120   return None;
0121 }
0122 
0123 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0124 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0125 
0126 //********************************************************************************//
0127 SiStripGainsCalibTreeWorker::SiStripGainsCalibTreeWorker(const edm::ParameterSet& iConfig) {
0128   MinTrackMomentum = iConfig.getUntrackedParameter<double>("minTrackMomentum", 3.0);
0129   MaxTrackMomentum = iConfig.getUntrackedParameter<double>("maxTrackMomentum", 99999.0);
0130   MinTrackEta = iConfig.getUntrackedParameter<double>("minTrackEta", -5.0);
0131   MaxTrackEta = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0);
0132   MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips", 2);
0133   MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits", 8);
0134   MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double>("MaxTrackChiOverNdf", 3);
0135   MaxTrackingIteration = iConfig.getUntrackedParameter<int>("MaxTrackingIteration", 7);
0136   AllowSaturation = iConfig.getUntrackedParameter<bool>("AllowSaturation", false);
0137   FirstSetOfConstants = iConfig.getUntrackedParameter<bool>("FirstSetOfConstants", true);
0138   Validation = iConfig.getUntrackedParameter<bool>("Validation", false);
0139   OldGainRemoving = iConfig.getUntrackedParameter<bool>("OldGainRemoving", false);
0140   useCalibration = iConfig.getUntrackedParameter<bool>("UseCalibration", false);
0141   doChargeMonitorPerPlane = iConfig.getUntrackedParameter<bool>("doChargeMonitorPerPlane", false);
0142   m_DQMdir = iConfig.getUntrackedParameter<std::string>("DQMdir", "AlCaReco/SiStripGains");
0143   m_calibrationMode = iConfig.getUntrackedParameter<std::string>("calibrationMode", "StdBunch");
0144   VChargeHisto = iConfig.getUntrackedParameter<std::vector<std::string>>("ChargeHisto");
0145 
0146   // fill in the mapping between the histogram indices and the (id,side,plane) tuple
0147   std::vector<std::pair<std::string, std::string>> hnames =
0148       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0149   for (unsigned int i = 0; i < hnames.size(); i++) {
0150     int id = APVGain::subdetectorId((hnames[i]).first);
0151     int side = APVGain::subdetectorSide((hnames[i]).first);
0152     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0153     int thick = APVGain::thickness((hnames[i]).first);
0154     std::string s = hnames[i].first;
0155 
0156     auto loc = APVloc(thick, id, side, plane, s);
0157     theTopologyMap.insert(std::make_pair(i, loc));
0158   }
0159 
0160   //Set the monitoring element tag and store
0161   dqm_tag_.reserve(7);
0162   dqm_tag_.clear();
0163   dqm_tag_.push_back("StdBunch");    // statistic collection from Standard Collision Bunch @ 3.8 T
0164   dqm_tag_.push_back("StdBunch0T");  // statistic collection from Standard Collision Bunch @ 0 T
0165   dqm_tag_.push_back("AagBunch");    // statistic collection from First Collision After Abort Gap @ 3.8 T
0166   dqm_tag_.push_back("AagBunch0T");  // statistic collection from First Collision After Abort Gap @ 0 T
0167   dqm_tag_.push_back("IsoMuon");     // statistic collection from Isolated Muon @ 3.8 T
0168   dqm_tag_.push_back("IsoMuon0T");   // statistic collection from Isolated Muon @ 0 T
0169   dqm_tag_.push_back("Harvest");     // statistic collection: Harvest
0170 
0171   hasProcessed_ = false;
0172 
0173   edm::ParameterSet swhallowgain_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("gain");
0174   CalibPrefix_ = swhallowgain_pset.getUntrackedParameter<std::string>("prefix");
0175   CalibSuffix_ = swhallowgain_pset.getUntrackedParameter<std::string>("suffix");
0176   edm::ParameterSet evtinfo_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("evtinfo");
0177   EventPrefix_ = evtinfo_pset.getUntrackedParameter<std::string>("prefix");
0178   EventSuffix_ = evtinfo_pset.getUntrackedParameter<std::string>("suffix");
0179   edm::ParameterSet track_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("tracks");
0180   TrackPrefix_ = track_pset.getUntrackedParameter<std::string>("prefix");
0181   TrackSuffix_ = track_pset.getUntrackedParameter<std::string>("suffix");
0182 
0183   calibTreeName_ = iConfig.getUntrackedParameter<std::string>("CalibTreeName");
0184   calibTreeFileNames_ = iConfig.getUntrackedParameter<std::vector<std::string>>("CalibTreeFiles");
0185 
0186   tTopoToken_ = esConsumes();
0187   tkGeomTokenBR_ = esConsumes<edm::Transition::BeginRun>();
0188   gainToken_ = esConsumes<edm::Transition::BeginRun>();
0189   qualityToken_ = esConsumes<edm::Transition::BeginRun>();
0190 }
0191 
0192 //********************************************************************************//
0193 void SiStripGainsCalibTreeWorker::dqmBeginRun(edm::Run const& run,
0194                                               edm::EventSetup const& iSetup,
0195                                               APVGain::APVGainHistograms& histograms) const {
0196   using namespace edm;
0197 
0198   // fills the APV collections at each begin run
0199   const TrackerGeometry* bareTkGeomPtr = &iSetup.getData(tkGeomTokenBR_);
0200   checkBookAPVColls(bareTkGeomPtr, histograms);
0201 
0202   const auto gainHandle = iSetup.getHandle(gainToken_);
0203   if (!gainHandle.isValid()) {
0204     edm::LogError("SiStripGainPCLWorker") << "gainHandle is not valid\n";
0205     exit(0);
0206   }
0207 
0208   const auto& SiStripQuality_ = &iSetup.getData(qualityToken_);
0209 
0210   for (unsigned int a = 0; a < histograms.APVsCollOrdered.size(); a++) {
0211     std::shared_ptr<stAPVGain> APV = histograms.APVsCollOrdered[a];
0212 
0213     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0214       continue;
0215 
0216     APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId, APV->APVId);
0217 
0218     if (gainHandle->getNumberOfTags() != 2) {
0219       edm::LogError("SiStripGainPCLWorker") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
0220       fflush(stdout);
0221       exit(0);
0222     };
0223     float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
0224     if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
0225       edm::LogWarning("SiStripGainPCLWorker") << "WARNING: ParticleGain in the global tag changed\n";
0226     APV->PreviousGain = newPreviousGain;
0227 
0228     float newPreviousGainTick = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
0229     if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
0230       edm::LogWarning("SiStripGainPCLWorker")
0231           << "WARNING: TickMarkGain in the global tag changed\n"
0232           << std::endl
0233           << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
0234           << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
0235           << std::endl;
0236     }
0237     APV->PreviousGainTick = newPreviousGainTick;
0238   }
0239 }
0240 
0241 //********************************************************************************//
0242 // ------------ method called for each event  ------------
0243 void SiStripGainsCalibTreeWorker::dqmAnalyze(edm::Event const& iEvent,
0244                                              edm::EventSetup const& iSetup,
0245                                              APVGain::APVGainHistograms const& histograms) const {
0246   if (!hasProcessed_) {
0247     const TrackerTopology* topo = &iSetup.getData(tTopoToken_);
0248 
0249     for (const auto& elem : theTopologyMap) {
0250       LogDebug("SiStripGainsCalibTreeWorker")
0251           << elem.first << " - " << elem.second.m_string << " " << elem.second.m_subdetectorId << " "
0252           << elem.second.m_subdetectorSide << " " << elem.second.m_subdetectorPlane << std::endl;
0253     }
0254 
0255     LogDebug("SiStripGainsCalibTreeWorker") << "for mode" << m_calibrationMode << std::endl;
0256 
0257     // data members
0258     unsigned int eventnumber = 0;
0259     unsigned int runnumber = 0;
0260 #ifdef ExtendedCALIBTree
0261     const std::vector<bool>* TrigTech = nullptr;
0262     const std::vector<double>* chargeoverpath = nullptr;
0263 #endif
0264     // Track data
0265     const std::vector<double>* trackchi2ndof = nullptr;
0266     const std::vector<float>* trackp = nullptr;
0267     const std::vector<float>* trackpt = nullptr;
0268     const std::vector<double>* tracketa = nullptr;
0269     const std::vector<double>* trackphi = nullptr;
0270     const std::vector<unsigned int>* trackhitsvalid = nullptr;
0271     const std::vector<int>* trackalgo = nullptr;
0272     // CalibTree data
0273     const std::vector<int>* trackindex = nullptr;
0274     const std::vector<unsigned int>* rawid = nullptr;
0275     const std::vector<double>* localdirx = nullptr;
0276     const std::vector<double>* localdiry = nullptr;
0277     const std::vector<double>* localdirz = nullptr;
0278     const std::vector<unsigned short>* firststrip = nullptr;
0279     const std::vector<unsigned short>* nstrips = nullptr;
0280     const std::vector<bool>* saturation = nullptr;
0281     const std::vector<bool>* overlapping = nullptr;
0282     const std::vector<bool>* farfromedge = nullptr;
0283     const std::vector<unsigned int>* charge = nullptr;
0284     const std::vector<double>* path = nullptr;
0285     const std::vector<unsigned char>* amplitude = nullptr;
0286     const std::vector<double>* gainused = nullptr;
0287     const std::vector<double>* gainusedTick = nullptr;
0288 
0289     TChain tree{calibTreeName_.c_str()};
0290     for (const auto& ctFn : calibTreeFileNames_) {
0291       tree.Add(ctFn.c_str());
0292     }
0293 
0294     tree.SetBranchAddress((EventPrefix_ + "event" + EventSuffix_).c_str(), &eventnumber, nullptr);
0295     tree.SetBranchAddress((EventPrefix_ + "run" + EventSuffix_).c_str(), &runnumber, nullptr);
0296 #ifdef ExtendedCALIBTree
0297     tree.SetBranchAddress((EventPrefix_ + "TrigTech" + EventSuffix_).c_str(), &TrigTech, nullptr);
0298     tree.SetBranchAddress((CalibPrefix_ + "chargeoverpath" + CalibSuffix_).c_str(), &chargeoverpath, nullptr);
0299 #endif
0300     tree.SetBranchAddress((TrackPrefix_ + "chi2ndof" + TrackSuffix_).c_str(), &trackchi2ndof, nullptr);
0301     tree.SetBranchAddress((TrackPrefix_ + "momentum" + TrackSuffix_).c_str(), &trackp, nullptr);
0302     tree.SetBranchAddress((TrackPrefix_ + "pt" + TrackSuffix_).c_str(), &trackpt, nullptr);
0303     tree.SetBranchAddress((TrackPrefix_ + "eta" + TrackSuffix_).c_str(), &tracketa, nullptr);
0304     tree.SetBranchAddress((TrackPrefix_ + "phi" + TrackSuffix_).c_str(), &trackphi, nullptr);
0305     tree.SetBranchAddress((TrackPrefix_ + "hitsvalid" + TrackSuffix_).c_str(), &trackhitsvalid, nullptr);
0306     tree.SetBranchAddress((TrackPrefix_ + "algo" + TrackSuffix_).c_str(), &trackalgo, nullptr);
0307 
0308     tree.SetBranchAddress((CalibPrefix_ + "trackindex" + CalibSuffix_).c_str(), &trackindex, nullptr);
0309     tree.SetBranchAddress((CalibPrefix_ + "rawid" + CalibSuffix_).c_str(), &rawid, nullptr);
0310     tree.SetBranchAddress((CalibPrefix_ + "localdirx" + CalibSuffix_).c_str(), &localdirx, nullptr);
0311     tree.SetBranchAddress((CalibPrefix_ + "localdiry" + CalibSuffix_).c_str(), &localdiry, nullptr);
0312     tree.SetBranchAddress((CalibPrefix_ + "localdirz" + CalibSuffix_).c_str(), &localdirz, nullptr);
0313     tree.SetBranchAddress((CalibPrefix_ + "firststrip" + CalibSuffix_).c_str(), &firststrip, nullptr);
0314     tree.SetBranchAddress((CalibPrefix_ + "nstrips" + CalibSuffix_).c_str(), &nstrips, nullptr);
0315     tree.SetBranchAddress((CalibPrefix_ + "saturation" + CalibSuffix_).c_str(), &saturation, nullptr);
0316     tree.SetBranchAddress((CalibPrefix_ + "overlapping" + CalibSuffix_).c_str(), &overlapping, nullptr);
0317     tree.SetBranchAddress((CalibPrefix_ + "farfromedge" + CalibSuffix_).c_str(), &farfromedge, nullptr);
0318     tree.SetBranchAddress((CalibPrefix_ + "charge" + CalibSuffix_).c_str(), &charge, nullptr);
0319     tree.SetBranchAddress((CalibPrefix_ + "path" + CalibSuffix_).c_str(), &path, nullptr);
0320 
0321     tree.SetBranchAddress((CalibPrefix_ + "amplitude" + CalibSuffix_).c_str(), &amplitude, nullptr);
0322     tree.SetBranchAddress((CalibPrefix_ + "gainused" + CalibSuffix_).c_str(), &gainused, nullptr);
0323     tree.SetBranchAddress((CalibPrefix_ + "gainusedTick" + CalibSuffix_).c_str(), &gainusedTick, nullptr);
0324 
0325     int elepos = statCollectionFromMode(m_calibrationMode.c_str());
0326     const auto nEntries = tree.GetEntries();
0327     printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0328     printf("Looping on the Tree          :");
0329     const auto treeProgStep = std::max(nEntries / 50, Long64_t(1));
0330     for (Long64_t i{0}; i != nEntries; ++i) {
0331       if ((i % treeProgStep) == 0) {
0332         printf(".");
0333         fflush(stdout);
0334       }
0335       tree.GetEntry(i);
0336 
0337       unsigned int FirstAmplitude = 0;
0338       for (unsigned int i = 0; i < charge->size(); i++) {
0339         FirstAmplitude += (*nstrips)[i];
0340         int TI = (*trackindex)[i];
0341 
0342         if ((*tracketa)[TI] < MinTrackEta)
0343           continue;
0344         if ((*tracketa)[TI] > MaxTrackEta)
0345           continue;
0346         if ((*trackp)[TI] < MinTrackMomentum)
0347           continue;
0348         if ((*trackp)[TI] > MaxTrackMomentum)
0349           continue;
0350         if ((*trackhitsvalid)[TI] < MinTrackHits)
0351           continue;
0352         if ((*trackchi2ndof)[TI] > MaxTrackChiOverNdf)
0353           continue;
0354         if ((*trackalgo)[TI] > MaxTrackingIteration)
0355           continue;
0356 
0357         std::shared_ptr<stAPVGain> APV = histograms.APVsColl.at(
0358             ((*rawid)[i] << 4) |
0359             ((*firststrip)[i] /
0360              128));  //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
0361 
0362         if (APV->SubDet > 2 && (*farfromedge)[i] == false)
0363           continue;
0364         if (APV->SubDet > 2 && (*overlapping)[i] == true)
0365           continue;
0366         if (APV->SubDet > 2 && (*saturation)[i] && !AllowSaturation)
0367           continue;
0368         if (APV->SubDet > 2 && (*nstrips)[i] > MaxNrStrips)
0369           continue;
0370 
0371         int Charge = 0;
0372         if (APV->SubDet > 2 && (useCalibration || !FirstSetOfConstants)) {
0373           bool Saturation = false;
0374           for (unsigned int s = 0; s < (*nstrips)[i]; s++) {
0375             int StripCharge = (*amplitude)[FirstAmplitude - (*nstrips)[i] + s];
0376             if (useCalibration && !FirstSetOfConstants) {
0377               StripCharge = (int)(StripCharge * (APV->PreviousGain / APV->CalibGain));
0378             } else if (useCalibration) {
0379               StripCharge = (int)(StripCharge / APV->CalibGain);
0380             } else if (!FirstSetOfConstants) {
0381               StripCharge = (int)(StripCharge * APV->PreviousGain);
0382             }
0383             if (StripCharge > 1024) {
0384               StripCharge = 255;
0385               Saturation = true;
0386             } else if (StripCharge > 254) {
0387               StripCharge = 254;
0388               Saturation = true;
0389             }
0390             Charge += StripCharge;
0391           }
0392           if (Saturation && !AllowSaturation)
0393             continue;
0394         } else if (APV->SubDet > 2) {
0395           Charge = (*charge)[i];
0396         } else {
0397           Charge = (*charge)[i] / 265.0;  //expected scale factor between pixel and strip charge
0398         }
0399 
0400         double ClusterChargeOverPath = ((double)Charge) / (*path)[i];
0401         if (APV->SubDet > 2) {
0402           if (Validation) {
0403             ClusterChargeOverPath /= (*gainused)[i];
0404           }
0405           if (OldGainRemoving) {
0406             ClusterChargeOverPath *= (*gainused)[i];
0407           }
0408         }
0409 
0410         // keep processing of pixel cluster charge until here
0411         if (APV->SubDet <= 2)
0412           continue;
0413 
0414         // real histogram for calibration
0415         histograms.Charge_Vs_Index[elepos]->Fill(APV->Index, ClusterChargeOverPath);
0416         LogDebug("SiStripGainsCalibTreeWorker")
0417             << " for mode " << m_calibrationMode << "\n"
0418             << " i " << i << " useCalibration " << useCalibration << " FirstSetOfConstants " << FirstSetOfConstants
0419             << " APV->PreviousGain " << APV->PreviousGain << " APV->CalibGain " << APV->CalibGain << " APV->DetId "
0420             << APV->DetId << " APV->Index " << APV->Index << " Charge " << Charge << " Path " << (*path)[i]
0421             << " ClusterChargeOverPath " << ClusterChargeOverPath << std::endl;
0422 
0423         // Fill monitoring histograms
0424         int mCharge1 = 0;
0425         int mCharge2 = 0;
0426         int mCharge3 = 0;
0427         int mCharge4 = 0;
0428         if (APV->SubDet > 2) {
0429           for (unsigned int s = 0; s < (*nstrips)[i]; s++) {
0430             int StripCharge = (*amplitude)[FirstAmplitude - (*nstrips)[i] + s];
0431             if (StripCharge > 1024)
0432               StripCharge = 255;
0433             else if (StripCharge > 254)
0434               StripCharge = 254;
0435             mCharge1 += StripCharge;
0436             mCharge2 += StripCharge;
0437             mCharge3 += StripCharge;
0438             mCharge4 += StripCharge;
0439           }
0440           // Revome gains for monitoring
0441           mCharge2 *= (*gainused)[i];                         // remove G2
0442           mCharge3 *= (*gainusedTick)[i];                     // remove G1
0443           mCharge4 *= ((*gainused)[i] * (*gainusedTick)[i]);  // remove G1 and G2
0444         }
0445 
0446         LogDebug("SiStripGainsCalibTreeWorker")
0447             << " full charge " << mCharge1 << " remove G2 " << mCharge2 << " remove G1 " << mCharge3 << " remove G1*G2 "
0448             << mCharge4 << std::endl;
0449 
0450         auto indices = APVGain::FetchIndices(theTopologyMap, (*rawid)[i], topo);
0451 
0452         for (auto m : indices)
0453           histograms.Charge_1[elepos][m]->Fill(((double)mCharge1) / (*path)[i]);
0454         for (auto m : indices)
0455           histograms.Charge_2[elepos][m]->Fill(((double)mCharge2) / (*path)[i]);
0456         for (auto m : indices)
0457           histograms.Charge_3[elepos][m]->Fill(((double)mCharge3) / (*path)[i]);
0458         for (auto m : indices)
0459           histograms.Charge_4[elepos][m]->Fill(((double)mCharge4) / (*path)[i]);
0460 
0461         if (APV->SubDet == StripSubdetector::TIB) {
0462           histograms.Charge_Vs_PathlengthTIB[elepos]->Fill((*path)[i], Charge);  // TIB
0463 
0464         } else if (APV->SubDet == StripSubdetector::TOB) {
0465           histograms.Charge_Vs_PathlengthTOB[elepos]->Fill((*path)[i], Charge);  // TOB
0466 
0467         } else if (APV->SubDet == StripSubdetector::TID) {
0468           if (APV->Eta < 0) {
0469             histograms.Charge_Vs_PathlengthTIDM[elepos]->Fill((*path)[i], Charge);
0470           }  // TID minus
0471           else if (APV->Eta > 0) {
0472             histograms.Charge_Vs_PathlengthTIDP[elepos]->Fill((*path)[i], Charge);
0473           }  // TID plus
0474 
0475         } else if (APV->SubDet == StripSubdetector::TEC) {
0476           if (APV->Eta < 0) {
0477             if (APV->Thickness < 0.04) {
0478               histograms.Charge_Vs_PathlengthTECM1[elepos]->Fill((*path)[i], Charge);
0479             }  // TEC minus, type 1
0480             else if (APV->Thickness > 0.04) {
0481               histograms.Charge_Vs_PathlengthTECM2[elepos]->Fill((*path)[i], Charge);
0482             }  // TEC minus, type 2
0483           } else if (APV->Eta > 0) {
0484             if (APV->Thickness < 0.04) {
0485               histograms.Charge_Vs_PathlengthTECP1[elepos]->Fill((*path)[i], Charge);
0486             }  // TEC plus, type 1
0487             else if (APV->Thickness > 0.04) {
0488               histograms.Charge_Vs_PathlengthTECP2[elepos]->Fill((*path)[i], Charge);
0489             }  // TEC plus, type 2
0490           }
0491         }
0492 
0493       }  // END OF ON-CLUSTER LOOP
0494 
0495       //LogDebug("SiStripGainsCalibTreeWorker")<<" for mode"<< m_calibrationMode
0496       //                   <<" entries in histogram:"<< histograms.Charge_Vs_Index[elepos].getEntries()
0497 
0498       //                   <<std::endl;
0499     }
0500 
0501     hasProcessed_ = true;
0502   }
0503 }
0504 
0505 //********************************************************************************//
0506 // ------------ method called once each job just before starting event loop  ------------
0507 void SiStripGainsCalibTreeWorker::checkBookAPVColls(const TrackerGeometry* bareTkGeomPtr,
0508                                                     APVGain::APVGainHistograms& histograms) const {
0509   if (bareTkGeomPtr) {  // pointer not yet set: called the first time => fill the APVColls
0510     auto const& Det = bareTkGeomPtr->dets();
0511 
0512     edm::LogInfo("SiStripGainsCalibTreeWorker") << " Resetting APV struct" << std::endl;
0513 
0514     unsigned int Index = 0;
0515 
0516     for (unsigned int i = 0; i < Det.size(); i++) {
0517       DetId Detid = Det[i]->geographicalId();
0518       int SubDet = Detid.subdetId();
0519 
0520       if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
0521           SubDet == StripSubdetector::TEC) {
0522         auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
0523         if (!DetUnit)
0524           continue;
0525 
0526         const StripTopology& Topo = DetUnit->specificTopology();
0527         unsigned int NAPV = Topo.nstrips() / 128;
0528 
0529         for (unsigned int j = 0; j < NAPV; j++) {
0530           auto APV = std::make_shared<stAPVGain>();
0531           APV->Index = Index;
0532           APV->Bin = -1;
0533           APV->DetId = Detid.rawId();
0534           APV->APVId = j;
0535           APV->SubDet = SubDet;
0536           APV->FitMPV = -1;
0537           APV->FitMPVErr = -1;
0538           APV->FitWidth = -1;
0539           APV->FitWidthErr = -1;
0540           APV->FitChi2 = -1;
0541           APV->FitNorm = -1;
0542           APV->Gain = -1;
0543           APV->PreviousGain = 1;
0544           APV->PreviousGainTick = 1;
0545           APV->x = DetUnit->position().basicVector().x();
0546           APV->y = DetUnit->position().basicVector().y();
0547           APV->z = DetUnit->position().basicVector().z();
0548           APV->Eta = DetUnit->position().basicVector().eta();
0549           APV->Phi = DetUnit->position().basicVector().phi();
0550           APV->R = DetUnit->position().basicVector().transverse();
0551           APV->Thickness = DetUnit->surface().bounds().thickness();
0552           APV->NEntries = 0;
0553           APV->isMasked = false;
0554 
0555           histograms.APVsCollOrdered.push_back(APV);
0556           histograms.APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0557           Index++;
0558           histograms.NStripAPVs++;
0559         }  // loop on APVs
0560       }    // if is Strips
0561     }      // loop on dets
0562 
0563     for (unsigned int i = 0; i < Det.size();
0564          i++) {  //Make two loop such that the Pixel information is added at the end --> make transition simpler
0565       DetId Detid = Det[i]->geographicalId();
0566       int SubDet = Detid.subdetId();
0567       if (SubDet == PixelSubdetector::PixelBarrel || SubDet == PixelSubdetector::PixelEndcap) {
0568         auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
0569         if (!DetUnit)
0570           continue;
0571 
0572         const PixelTopology& Topo = DetUnit->specificTopology();
0573         unsigned int NROCRow = Topo.nrows() / (80.);
0574         unsigned int NROCCol = Topo.ncolumns() / (52.);
0575 
0576         for (unsigned int j = 0; j < NROCRow; j++) {
0577           for (unsigned int i = 0; i < NROCCol; i++) {
0578             auto APV = std::make_shared<stAPVGain>();
0579             APV->Index = Index;
0580             APV->Bin = -1;
0581             APV->DetId = Detid.rawId();
0582             APV->APVId = (j << 3 | i);
0583             APV->SubDet = SubDet;
0584             APV->FitMPV = -1;
0585             APV->FitMPVErr = -1;
0586             APV->FitWidth = -1;
0587             APV->FitWidthErr = -1;
0588             APV->FitChi2 = -1;
0589             APV->Gain = -1;
0590             APV->PreviousGain = 1;
0591             APV->PreviousGainTick = 1;
0592             APV->x = DetUnit->position().basicVector().x();
0593             APV->y = DetUnit->position().basicVector().y();
0594             APV->z = DetUnit->position().basicVector().z();
0595             APV->Eta = DetUnit->position().basicVector().eta();
0596             APV->Phi = DetUnit->position().basicVector().phi();
0597             APV->R = DetUnit->position().basicVector().transverse();
0598             APV->Thickness = DetUnit->surface().bounds().thickness();
0599             APV->isMasked = false;  //SiPixelQuality_->IsModuleBad(Detid.rawId());
0600             APV->NEntries = 0;
0601 
0602             histograms.APVsCollOrdered.push_back(APV);
0603             histograms.APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0604             Index++;
0605             histograms.NPixelDets++;
0606 
0607           }  // loop on ROC cols
0608         }    // loop on ROC rows
0609       }      // if Pixel
0610     }        // loop on Dets
0611   }          //if (!bareTkGeomPtr_) ...
0612 }
0613 
0614 //********************************************************************************//
0615 void SiStripGainsCalibTreeWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0616   edm::ParameterSetDescription desc;
0617   desc.setUnknown();
0618   descriptions.addDefault(desc);
0619 }
0620 
0621 //********************************************************************************//
0622 void SiStripGainsCalibTreeWorker::bookHistograms(DQMStore::IBooker& ibooker,
0623                                                  edm::Run const& run,
0624                                                  edm::EventSetup const& setup,
0625                                                  APVGain::APVGainHistograms& histograms) const {
0626   ibooker.cd();
0627   std::string dqm_dir = m_DQMdir;
0628   const char* tag = dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str();
0629 
0630   edm::LogInfo("SiStripGainsCalibTreeWorker")
0631       << "Setting " << dqm_dir << " in DQM and booking histograms for tag " << tag << std::endl;
0632 
0633   ibooker.setCurrentFolder(dqm_dir);
0634 
0635   std::string stag(tag);
0636   if (!stag.empty() && stag[0] != '_')
0637     stag.insert(0, 1, '_');
0638 
0639   std::string cvi = std::string("Charge_Vs_Index") + stag;
0640   std::string cvpTIB = std::string("Charge_Vs_PathlengthTIB") + stag;
0641   std::string cvpTOB = std::string("Charge_Vs_PathlengthTOB") + stag;
0642   std::string cvpTIDP = std::string("Charge_Vs_PathlengthTIDP") + stag;
0643   std::string cvpTIDM = std::string("Charge_Vs_PathlengthTIDM") + stag;
0644   std::string cvpTECP1 = std::string("Charge_Vs_PathlengthTECP1") + stag;
0645   std::string cvpTECP2 = std::string("Charge_Vs_PathlengthTECP2") + stag;
0646   std::string cvpTECM1 = std::string("Charge_Vs_PathlengthTECM1") + stag;
0647   std::string cvpTECM2 = std::string("Charge_Vs_PathlengthTECM2") + stag;
0648 
0649   int elepos = statCollectionFromMode(tag);
0650 
0651   histograms.Charge_Vs_Index.reserve(dqm_tag_.size());
0652   histograms.Charge_Vs_PathlengthTIB.reserve(dqm_tag_.size());
0653   histograms.Charge_Vs_PathlengthTOB.reserve(dqm_tag_.size());
0654   histograms.Charge_Vs_PathlengthTIDP.reserve(dqm_tag_.size());
0655   histograms.Charge_Vs_PathlengthTIDM.reserve(dqm_tag_.size());
0656   histograms.Charge_Vs_PathlengthTECP1.reserve(dqm_tag_.size());
0657   histograms.Charge_Vs_PathlengthTECP2.reserve(dqm_tag_.size());
0658   histograms.Charge_Vs_PathlengthTECM1.reserve(dqm_tag_.size());
0659   histograms.Charge_Vs_PathlengthTECM2.reserve(dqm_tag_.size());
0660 
0661   // The cluster charge is stored by exploiting a non uniform binning in order
0662   // reduce the histogram memory size. The bin width is relaxed with a falling
0663   // exponential function and the bin boundaries are stored in the binYarray.
0664   // The binXarray is used to provide as many bins as the APVs.
0665   //
0666   // More details about this implementations are here:
0667   // https://indico.cern.ch/event/649344/contributions/2672267/attachments/1498323/2332518/OptimizeChHisto.pdf
0668 
0669   std::vector<float> binXarray;
0670   binXarray.reserve(histograms.NStripAPVs + 1);
0671   for (unsigned int a = 0; a <= histograms.NStripAPVs; a++) {
0672     binXarray.push_back((float)a);
0673   }
0674 
0675   std::array<float, 688> binYarray;
0676   double p0 = 5.445;
0677   double p1 = 0.002113;
0678   double p2 = 69.01576;
0679   double y = 0.;
0680   for (int b = 0; b < 687; b++) {
0681     binYarray[b] = y;
0682     if (y <= 902.)
0683       y = y + 2.;
0684     else
0685       y = (p0 - log(exp(p0 - p1 * y) - p2 * p1)) / p1;
0686   }
0687   binYarray[687] = 4000.;
0688 
0689   histograms.Charge_1[elepos].clear();
0690   histograms.Charge_2[elepos].clear();
0691   histograms.Charge_3[elepos].clear();
0692   histograms.Charge_4[elepos].clear();
0693 
0694   auto it = histograms.Charge_Vs_Index.begin();
0695   histograms.Charge_Vs_Index.insert(
0696       it + elepos,
0697       ibooker.book2S(cvi.c_str(), cvi.c_str(), histograms.NStripAPVs, &binXarray[0], 687, binYarray.data()));
0698 
0699   it = histograms.Charge_Vs_PathlengthTIB.begin();
0700   histograms.Charge_Vs_PathlengthTIB.insert(it + elepos,
0701                                             ibooker.book2S(cvpTIB.c_str(), cvpTIB.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0702 
0703   it = histograms.Charge_Vs_PathlengthTOB.begin();
0704   histograms.Charge_Vs_PathlengthTOB.insert(it + elepos,
0705                                             ibooker.book2S(cvpTOB.c_str(), cvpTOB.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0706 
0707   it = histograms.Charge_Vs_PathlengthTIDP.begin();
0708   histograms.Charge_Vs_PathlengthTIDP.insert(
0709       it + elepos, ibooker.book2S(cvpTIDP.c_str(), cvpTIDP.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0710 
0711   it = histograms.Charge_Vs_PathlengthTIDM.begin();
0712   histograms.Charge_Vs_PathlengthTIDM.insert(
0713       it + elepos, ibooker.book2S(cvpTIDM.c_str(), cvpTIDM.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0714 
0715   it = histograms.Charge_Vs_PathlengthTECP1.begin();
0716   histograms.Charge_Vs_PathlengthTECP1.insert(
0717       it + elepos, ibooker.book2S(cvpTECP1.c_str(), cvpTECP1.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0718 
0719   it = histograms.Charge_Vs_PathlengthTECP2.begin();
0720   histograms.Charge_Vs_PathlengthTECP2.insert(
0721       it + elepos, ibooker.book2S(cvpTECP2.c_str(), cvpTECP2.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0722 
0723   it = histograms.Charge_Vs_PathlengthTECM1.begin();
0724   histograms.Charge_Vs_PathlengthTECM1.insert(
0725       it + elepos, ibooker.book2S(cvpTECM1.c_str(), cvpTECM1.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0726 
0727   it = histograms.Charge_Vs_PathlengthTECM2.begin();
0728   histograms.Charge_Vs_PathlengthTECM2.insert(
0729       it + elepos, ibooker.book2S(cvpTECM2.c_str(), cvpTECM2.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0730 
0731   std::vector<std::pair<std::string, std::string>> hnames =
0732       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0733   for (unsigned int i = 0; i < hnames.size(); i++) {
0734     std::string htag = (hnames[i]).first + stag;
0735     histograms.Charge_1[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0736   }
0737 
0738   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG2");
0739   for (unsigned int i = 0; i < hnames.size(); i++) {
0740     std::string htag = (hnames[i]).first + stag;
0741     histograms.Charge_2[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0742   }
0743 
0744   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1");
0745   for (unsigned int i = 0; i < hnames.size(); i++) {
0746     std::string htag = (hnames[i]).first + stag;
0747     histograms.Charge_3[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0748   }
0749 
0750   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1G2");
0751   for (unsigned int i = 0; i < hnames.size(); i++) {
0752     std::string htag = (hnames[i]).first + stag;
0753     histograms.Charge_4[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0754   }
0755 }
0756 
0757 #include "FWCore/PluginManager/interface/ModuleDef.h"
0758 DEFINE_FWK_MODULE(SiStripGainsCalibTreeWorker);