Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:59:40

0001 // Original Author:  Loic QUERTENMONT
0002 //         Created:  Mon Nov  16 08:55:18 CET 2009
0003 
0004 #include <memory>
0005 #include <iostream>
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 #include "FWCore/Utilities/interface/EDGetToken.h"
0014 
0015 #include "CondFormats/RunInfo/interface/RunInfo.h"
0016 
0017 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0018 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0019 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0020 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0021 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0022 
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0025 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0026 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0027 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0028 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0029 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0030 
0031 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
0032 
0033 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0034 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
0035 
0036 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
0037 #include "DataFormats/TrackReco/interface/Track.h"
0038 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0039 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0040 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0041 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0042 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0043 #include "DataFormats/DetId/interface/DetId.h"
0044 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0045 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0046 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
0047 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0048 
0049 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
0050 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
0051 
0052 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0053 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0054 
0055 #include "DQMServices/Core/interface/DQMStore.h"
0056 #include "FWCore/ServiceRegistry/interface/Service.h"
0057 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0058 
0059 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0060 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
0061 
0062 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0063 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0064 
0065 #include "TFile.h"
0066 #include "TObjString.h"
0067 #include "TString.h"
0068 #include "TH1F.h"
0069 #include "TH2S.h"
0070 #include "TProfile.h"
0071 #include "TF1.h"
0072 #include "TROOT.h"
0073 #include "TTree.h"
0074 #include "TChain.h"
0075 
0076 // user includes
0077 #include "CalibTracker/SiStripChannelGain/interface/APVGainStruct.h"
0078 #include "CalibTracker/SiStripChannelGain/interface/APVGainHelpers.h"
0079 
0080 #include <unordered_map>
0081 #include <array>
0082 
0083 using namespace edm;
0084 using namespace reco;
0085 using namespace std;
0086 using namespace APVGain;
0087 
0088 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
0089 public:
0090   typedef dqm::legacy::MonitorElement MonitorElement;
0091   typedef dqm::legacy::DQMStore DQMStore;
0092   explicit SiStripGainFromCalibTree(const edm::ParameterSet&);
0093   ~SiStripGainFromCalibTree() override;
0094 
0095 private:
0096   void algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
0097   void algoEndRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
0098   void algoBeginJob(const edm::EventSetup& iSetup) override;
0099   void algoEndJob() override;
0100   void algoAnalyze(const edm::Event&, const edm::EventSetup&) override;
0101 
0102   int statCollectionFromMode(const char* tag) const;
0103   void bookDQMHistos(const char* dqm_dir, const char* tag);
0104 
0105   bool isBFieldConsistentWithMode(const edm::EventSetup& iSetup) const;
0106   void swapBFieldMode(void);
0107 
0108   void merge(TH2* A, TH2* B);  //needed to add histograms with different number of bins
0109   void algoAnalyzeTheTree();
0110   void algoComputeMPVandGain();
0111   void processEvent();  //what really does the job
0112 
0113   void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange = 50, double HighRange = 5400);
0114   bool IsGoodLandauFit(double* FitResults);
0115   void storeOnTree(TFileService* tfs);
0116   void qualityMonitor();
0117   void MakeCalibrationMap();
0118   bool produceTagFilter();
0119 
0120   template <typename T>
0121   inline edm::Handle<T> connect(const T*& ptr, edm::EDGetTokenT<T> token, const edm::Event& evt) {
0122     edm::Handle<T> handle;
0123     evt.getByToken(token, handle);
0124     ptr = handle.product();
0125     return handle;  //return handle to keep alive pointer (safety first)
0126   }
0127 
0128   std::unique_ptr<SiStripApvGain> getNewObject() override;
0129 
0130   TFileService* tfs;
0131   DQMStore* dbe;
0132   double MagFieldCurrentTh;
0133   double MinNrEntries;
0134   double MaxMPVError;
0135   double MaxChi2OverNDF;
0136   double MinTrackMomentum;
0137   double MaxTrackMomentum;
0138   double MinTrackEta;
0139   double MaxTrackEta;
0140   unsigned int MaxNrStrips;
0141   unsigned int MinTrackHits;
0142   double MaxTrackChiOverNdf;
0143   int MaxTrackingIteration;
0144   bool AllowSaturation;
0145   bool FirstSetOfConstants;
0146   bool Validation;
0147   bool OldGainRemoving;
0148   int CalibrationLevel;
0149 
0150   bool saveSummary;
0151   bool useCalibration;
0152   bool m_harvestingMode;
0153   bool m_splitDQMstat;
0154   bool doChargeMonitorPerPlane;  /*!< Charge monitor per detector plane */
0155   std::string m_calibrationMode; /*!< Type of statistics for the calibration */
0156   std::string m_calibrationPath;
0157   std::string m_DQMdir;                  /*!< DQM folder hosting the charge statistics and the monitor plots */
0158   std::vector<std::string> VChargeHisto; /*!< Charge monitor plots to be output */
0159 
0160   double tagCondition_NClusters;
0161   double tagCondition_GoodFrac;
0162 
0163   std::string AlgoMode;
0164   std::string OutputGains;
0165   vector<string> VInputFiles;
0166 
0167   //enum statistic_type {None=-1, StdBunch, StdBunch0T, FaABunch, FaABunch0T, IsoBunch, IsoBunch0T, Harvest};
0168 
0169   std::vector<string> dqm_tag_;
0170   std::string booked_dir_;
0171 
0172   std::vector<MonitorElement*> Charge_Vs_Index;         /*!< Charge per cm for each detector id */
0173   std::array<std::vector<APVGain::APVmon>, 7> Charge_1; /*!< Charge per cm per layer / wheel */
0174   std::array<std::vector<APVGain::APVmon>, 7> Charge_2; /*!< Charge per cm per layer / wheel without G2 */
0175   std::array<std::vector<APVGain::APVmon>, 7> Charge_3; /*!< Charge per cm per layer / wheel without G1 */
0176   std::array<std::vector<APVGain::APVmon>, 7> Charge_4; /*!< Charge per cm per layer / wheel without G1 and G1*/
0177 
0178   std::vector<MonitorElement*> Charge_Vs_PathlengthTIB;   /*!< Charge vs pathlength in TIB */
0179   std::vector<MonitorElement*> Charge_Vs_PathlengthTOB;   /*!< Charge vs pathlength in TOB */
0180   std::vector<MonitorElement*> Charge_Vs_PathlengthTIDP;  /*!< Charge vs pathlength in TIDP */
0181   std::vector<MonitorElement*> Charge_Vs_PathlengthTIDM;  /*!< Charge vs pathlength in TIDM */
0182   std::vector<MonitorElement*> Charge_Vs_PathlengthTECP1; /*!< Charge vs pathlength in TECP thin */
0183   std::vector<MonitorElement*> Charge_Vs_PathlengthTECP2; /*!< Charge vs pathlength in TECP thick */
0184   std::vector<MonitorElement*> Charge_Vs_PathlengthTECM1; /*!< Charge vs pathlength in TECP thin */
0185   std::vector<MonitorElement*> Charge_Vs_PathlengthTECM2; /*!< Charge vs pathlength in TECP thick */
0186 
0187   //std::vector<MonitorElement*>  Charge_Vs_Index_Absolute;
0188 
0189   //Validation histograms
0190   MonitorElement* MPV_Vs_EtaTIB;      /*!< MPV vs Eta for TIB planes */
0191   MonitorElement* MPV_Vs_EtaTID;      /*!< MPV vs Eta for TID planes */
0192   MonitorElement* MPV_Vs_EtaTOB;      /*!< MPV vs Eta for TOB planes */
0193   MonitorElement* MPV_Vs_EtaTEC;      /*!< MPV vs Eta for TEC planes */
0194   MonitorElement* MPV_Vs_EtaTECthin;  /*!< MPV vs Eta for TEC thin planes */
0195   MonitorElement* MPV_Vs_EtaTECthick; /*!< MPV vs Eta for TEC tick planes */
0196 
0197   MonitorElement* MPV_Vs_PhiTIB;      /*!< MPV vs Phi for TIB planes */
0198   MonitorElement* MPV_Vs_PhiTID;      /*!< MPV vs Phi for TID planes */
0199   MonitorElement* MPV_Vs_PhiTOB;      /*!< MPV vs Phi for TOB planes */
0200   MonitorElement* MPV_Vs_PhiTEC;      /*!< MPV vs Phi for TEC planes */
0201   MonitorElement* MPV_Vs_PhiTECthin;  /*!< MPV vs Phi for TEC thin planes */
0202   MonitorElement* MPV_Vs_PhiTECthick; /*!< MPV vs Phi for TID tick planes */
0203 
0204   MonitorElement* NoMPVmasked; /*!< R,Z map of missing APV calibration (masked modules)*/
0205   MonitorElement* NoMPVfit;    /*!< R,Z map of missing APV calibration (no fit) */
0206 
0207   MonitorElement* Gains;        /*!< distribution of gain factors */
0208   MonitorElement* MPVs;         /*!< distribution of MPVs */
0209   MonitorElement* MPVs320;      /*!< distribution of MPVs for thin sensors */
0210   MonitorElement* MPVs500;      /*!< distribution of MPVs for tick sensors */
0211   MonitorElement* MPVsTIB;      /*!< distribution of MPVs for TIB planes */
0212   MonitorElement* MPVsTID;      /*!< distribution of MPVs for TID disks */
0213   MonitorElement* MPVsTIDP;     /*!< distribution of MPVs for TIDP disks */
0214   MonitorElement* MPVsTIDM;     /*!< distribution of MPVs for TIDM disks */
0215   MonitorElement* MPVsTOB;      /*!< distribution of MPVs for TOB planes */
0216   MonitorElement* MPVsTEC;      /*!< distribution of MPVs for TEC disks */
0217   MonitorElement* MPVsTECP;     /*!< distribution of MPVs for TECP disks */
0218   MonitorElement* MPVsTECM;     /*!< distribution of MPVs for TECM disks */
0219   MonitorElement* MPVsTECthin;  /*!< distribution of MPVs for TEC thin sensors */
0220   MonitorElement* MPVsTECthick; /*!< distribution of MPVs for TEC tick sensors */
0221   MonitorElement* MPVsTECP1;    /*!< distribution of MPVs for TECP thin sensors */
0222   MonitorElement* MPVsTECP2;    /*!< distribution of MPVs for TECP tick sensors */
0223   MonitorElement* MPVsTECM1;    /*!< distribution of MPVs for TECM thin sensors */
0224   MonitorElement* MPVsTECM2;    /*!< distribution of MPVs for TECM tick sensors */
0225 
0226   MonitorElement* MPVError;      /*!< error of Landau fit */
0227   MonitorElement* MPVErrorVsMPV; /*!< error of Landau fit vs MPV */
0228   MonitorElement* MPVErrorVsEta; /*!< error of Landau fit vs Eta */
0229   MonitorElement* MPVErrorVsPhi; /*!< error of Landau fit vs Phi */
0230   MonitorElement* MPVErrorVsN;   /*!< error of landau fit vs number of entries */
0231 
0232   MonitorElement* DiffWRTPrevGainTIB; /*!< ratio Gain / PreviousGain for TIB layers */
0233   MonitorElement* DiffWRTPrevGainTID; /*!< ratio Gain / PreviousGain for TID disks */
0234   MonitorElement* DiffWRTPrevGainTOB; /*!< ratio Gain / PreviousGain for TOB layers */
0235   MonitorElement* DiffWRTPrevGainTEC; /*!< ratio gain / PreviousGain for TEC disks */
0236 
0237   MonitorElement* GainVsPrevGainTIB; /*!< Gain vs PreviousGain for TIB */
0238   MonitorElement* GainVsPrevGainTID; /*!< Gain vs PreviousGain for TID */
0239   MonitorElement* GainVsPrevGainTOB; /*!< Gain vs PreviousGain for TOB */
0240   MonitorElement* GainVsPrevGainTEC; /*!< Gain vs PreviousGain for TEC */
0241 
0242   std::vector<APVGain::APVmon> newCharge;
0243 
0244   unsigned int NEvent;
0245   unsigned int NTrack;
0246   unsigned int NClusterStrip;
0247   unsigned int NClusterPixel;
0248   int NStripAPVs;
0249   int NPixelDets;
0250   unsigned int SRun;
0251   unsigned int ERun;
0252   unsigned int GOOD;
0253   unsigned int BAD;
0254   unsigned int MASKED;
0255 
0256   //Data members for processing
0257 
0258   //Event data
0259   unsigned int eventnumber = 0;
0260   unsigned int runnumber = 0;
0261   const std::vector<bool>* TrigTech = nullptr;
0262   edm::EDGetTokenT<std::vector<bool>> TrigTech_token_;
0263 
0264   // Track data
0265   const std::vector<double>* trackchi2ndof = nullptr;
0266   edm::EDGetTokenT<std::vector<double>> trackchi2ndof_token_;
0267   const std::vector<float>* trackp = nullptr;
0268   edm::EDGetTokenT<std::vector<float>> trackp_token_;
0269   const std::vector<float>* trackpt = nullptr;
0270   edm::EDGetTokenT<std::vector<float>> trackpt_token_;
0271   const std::vector<double>* tracketa = nullptr;
0272   edm::EDGetTokenT<std::vector<double>> tracketa_token_;
0273   const std::vector<double>* trackphi = nullptr;
0274   edm::EDGetTokenT<std::vector<double>> trackphi_token_;
0275   const std::vector<unsigned int>* trackhitsvalid = nullptr;
0276   edm::EDGetTokenT<std::vector<unsigned int>> trackhitsvalid_token_;
0277   const std::vector<int>* trackalgo = nullptr;
0278   edm::EDGetTokenT<std::vector<int>> trackalgo_token_;
0279 
0280   // CalibTree data
0281   const std::vector<int>* trackindex = nullptr;
0282   edm::EDGetTokenT<std::vector<int>> trackindex_token_;
0283   const std::vector<unsigned int>* rawid = nullptr;
0284   edm::EDGetTokenT<std::vector<unsigned int>> rawid_token_;
0285   const std::vector<double>* localdirx = nullptr;
0286   edm::EDGetTokenT<std::vector<double>> localdirx_token_;
0287   const std::vector<double>* localdiry = nullptr;
0288   edm::EDGetTokenT<std::vector<double>> localdiry_token_;
0289   const std::vector<double>* localdirz = nullptr;
0290   edm::EDGetTokenT<std::vector<double>> localdirz_token_;
0291   const std::vector<unsigned short>* firststrip = nullptr;
0292   edm::EDGetTokenT<std::vector<unsigned short>> firststrip_token_;
0293   const std::vector<unsigned short>* nstrips = nullptr;
0294   edm::EDGetTokenT<std::vector<unsigned short>> nstrips_token_;
0295   const std::vector<bool>* saturation = nullptr;
0296   edm::EDGetTokenT<std::vector<bool>> saturation_token_;
0297   const std::vector<bool>* overlapping = nullptr;
0298   edm::EDGetTokenT<std::vector<bool>> overlapping_token_;
0299   const std::vector<bool>* farfromedge = nullptr;
0300   edm::EDGetTokenT<std::vector<bool>> farfromedge_token_;
0301   const std::vector<unsigned int>* charge = nullptr;
0302   edm::EDGetTokenT<std::vector<unsigned int>> charge_token_;
0303   const std::vector<double>* path = nullptr;
0304   edm::EDGetTokenT<std::vector<double>> path_token_;
0305   const std::vector<double>* chargeoverpath = nullptr;
0306   edm::EDGetTokenT<std::vector<double>> chargeoverpath_token_;
0307   const std::vector<unsigned char>* amplitude = nullptr;
0308   edm::EDGetTokenT<std::vector<unsigned char>> amplitude_token_;
0309   const std::vector<double>* gainused = nullptr;
0310   edm::EDGetTokenT<std::vector<double>> gainused_token_;
0311   const std::vector<double>* gainusedTick = nullptr;
0312   edm::EDGetTokenT<std::vector<double>> gainusedTick_token_;
0313 
0314   string EventPrefix_;  //("");
0315   string EventSuffix_;  //("");
0316   string TrackPrefix_;  //("track");
0317   string TrackSuffix_;  //("");
0318   string CalibPrefix_;  //("GainCalibration");
0319   string CalibSuffix_;  //("");
0320 
0321 private:
0322   std::vector<stAPVGain*> APVsCollOrdered;
0323   std::unordered_map<unsigned int, stAPVGain*> APVsColl;
0324 
0325   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0326   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0327   edm::ESGetToken<RunInfo, RunInfoRcd> runInfoToken_;
0328   edm::ESGetToken<SiStripGain, SiStripGainRcd> gainToken_;
0329   edm::ESGetToken<SiStripQuality, SiStripQualityRcd> qualityToken_;
0330   const TrackerTopology* tTopo_ = nullptr;
0331 };
0332 
0333 inline int SiStripGainFromCalibTree::statCollectionFromMode(const char* tag) const {
0334   std::vector<string>::const_iterator it = dqm_tag_.begin();
0335   while (it != dqm_tag_.end()) {
0336     if (*it == std::string(tag))
0337       return it - dqm_tag_.begin();
0338     it++;
0339   }
0340 
0341   if (std::string(tag).empty())
0342     return 0;  // return StdBunch calibration mode for backward compatibility
0343 
0344   return None;
0345 }
0346 
0347 void SiStripGainFromCalibTree::merge(TH2* A, TH2* B) {
0348   if (A->GetNbinsX() == B->GetNbinsX()) {
0349     A->Add(B);
0350   } else {
0351     for (int x = 0; x <= B->GetNbinsX() + 1; x++) {
0352       for (int y = 0; y <= B->GetNbinsY() + 1; y++) {
0353         A->SetBinContent(x, y, A->GetBinContent(x, y) + B->GetBinContent(x, y));
0354       }
0355     }
0356   }
0357 }
0358 
0359 SiStripGainFromCalibTree::SiStripGainFromCalibTree(const edm::ParameterSet& iConfig)
0360     : ConditionDBWriter<SiStripApvGain>(iConfig) {
0361   usesResource(TFileService::kSharedResource);
0362   OutputGains = iConfig.getParameter<std::string>("OutputGains");
0363   AlgoMode = iConfig.getUntrackedParameter<std::string>("AlgoMode", "CalibTree");
0364   MagFieldCurrentTh = iConfig.getUntrackedParameter<double>("MagFieldCurrentTh", 2000.);
0365   MinNrEntries = iConfig.getUntrackedParameter<double>("minNrEntries", 20);
0366   MaxMPVError = iConfig.getUntrackedParameter<double>("maxMPVError", 500.0);
0367   MaxChi2OverNDF = iConfig.getUntrackedParameter<double>("maxChi2OverNDF", 5.0);
0368   MinTrackMomentum = iConfig.getUntrackedParameter<double>("minTrackMomentum", 3.0);
0369   MaxTrackMomentum = iConfig.getUntrackedParameter<double>("maxTrackMomentum", 99999.0);
0370   MinTrackEta = iConfig.getUntrackedParameter<double>("minTrackEta", -5.0);
0371   MaxTrackEta = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0);
0372   MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips", 2);
0373   MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits", 8);
0374   MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double>("MaxTrackChiOverNdf", 3);
0375   MaxTrackingIteration = iConfig.getUntrackedParameter<int>("MaxTrackingIteration", 7);
0376   AllowSaturation = iConfig.getUntrackedParameter<bool>("AllowSaturation", false);
0377   FirstSetOfConstants = iConfig.getUntrackedParameter<bool>("FirstSetOfConstants", true);
0378   Validation = iConfig.getUntrackedParameter<bool>("Validation", false);
0379   OldGainRemoving = iConfig.getUntrackedParameter<bool>("OldGainRemoving", false);
0380 
0381   CalibrationLevel = iConfig.getUntrackedParameter<int>("CalibrationLevel", 0);
0382   VInputFiles = iConfig.getUntrackedParameter<vector<string>>("InputFiles");
0383   VChargeHisto = iConfig.getUntrackedParameter<vector<string>>("ChargeHisto");
0384 
0385   useCalibration = iConfig.getUntrackedParameter<bool>("UseCalibration", false);
0386   m_harvestingMode = iConfig.getUntrackedParameter<bool>("harvestingMode", false);
0387   m_splitDQMstat = iConfig.getUntrackedParameter<bool>("splitDQMstat", false);
0388   m_calibrationMode = iConfig.getUntrackedParameter<string>("calibrationMode", "StdBunch");
0389   m_calibrationPath = iConfig.getUntrackedParameter<string>("calibrationPath");
0390   m_DQMdir = iConfig.getUntrackedParameter<string>("DQMdir", "AlCaReco/SiStripGains");
0391 
0392   tagCondition_NClusters = iConfig.getUntrackedParameter<double>("NClustersForTagProd", 2E8);
0393   tagCondition_GoodFrac = iConfig.getUntrackedParameter<double>("GoodFracForTagProd", 0.95);
0394 
0395   saveSummary = iConfig.getUntrackedParameter<bool>("saveSummary", false);
0396 
0397   doChargeMonitorPerPlane = iConfig.getUntrackedParameter<bool>("doChargeMonitorPerPlane", false);
0398 
0399   // Gather DQM Service
0400   dbe = edm::Service<DQMStore>().operator->();
0401 
0402   //Set the monitoring element tag and store
0403   dqm_tag_.reserve(7);
0404   dqm_tag_.clear();
0405   dqm_tag_.push_back("StdBunch");    // statistic collection from Standard Collision Bunch @ 3.8 T
0406   dqm_tag_.push_back("StdBunch0T");  // statistic collection from Standard Collision Bunch @ 0 T
0407   dqm_tag_.push_back("AagBunch");    // statistic collection from First Collision After Abort Gap @ 3.8 T
0408   dqm_tag_.push_back("AagBunch0T");  // statistic collection from First Collision After Abort Gap @ 0 T
0409   dqm_tag_.push_back("IsoMuon");     // statistic collection from Isolated Muon @ 3.8 T
0410   dqm_tag_.push_back("IsoMuon0T");   // statistic collection from Isolated Muon @ 0 T
0411   dqm_tag_.push_back("Harvest");     // statistic collection: Harvest
0412 
0413   Charge_Vs_Index.insert(Charge_Vs_Index.begin(), dqm_tag_.size(), nullptr);
0414   //Charge_Vs_Index_Absolute.insert( Charge_Vs_Index_Absolute.begin(), dqm_tag_.size(), 0);
0415   Charge_Vs_PathlengthTIB.insert(Charge_Vs_PathlengthTIB.begin(), dqm_tag_.size(), nullptr);
0416   Charge_Vs_PathlengthTOB.insert(Charge_Vs_PathlengthTOB.begin(), dqm_tag_.size(), nullptr);
0417   Charge_Vs_PathlengthTIDP.insert(Charge_Vs_PathlengthTIDP.begin(), dqm_tag_.size(), nullptr);
0418   Charge_Vs_PathlengthTIDM.insert(Charge_Vs_PathlengthTIDM.begin(), dqm_tag_.size(), nullptr);
0419   Charge_Vs_PathlengthTECP1.insert(Charge_Vs_PathlengthTECP1.begin(), dqm_tag_.size(), nullptr);
0420   Charge_Vs_PathlengthTECP2.insert(Charge_Vs_PathlengthTECP2.begin(), dqm_tag_.size(), nullptr);
0421   Charge_Vs_PathlengthTECM1.insert(Charge_Vs_PathlengthTECM1.begin(), dqm_tag_.size(), nullptr);
0422   Charge_Vs_PathlengthTECM2.insert(Charge_Vs_PathlengthTECM2.begin(), dqm_tag_.size(), nullptr);
0423 
0424   // configure token for gathering the ntuple variables
0425   edm::ParameterSet swhallowgain_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("gain");
0426 
0427   string label = swhallowgain_pset.getUntrackedParameter<string>("label");
0428   CalibPrefix_ = swhallowgain_pset.getUntrackedParameter<string>("prefix");
0429   CalibSuffix_ = swhallowgain_pset.getUntrackedParameter<string>("suffix");
0430 
0431   trackindex_token_ = consumes<std::vector<int>>(edm::InputTag(label, CalibPrefix_ + "trackindex" + CalibSuffix_));
0432   rawid_token_ = consumes<std::vector<unsigned int>>(edm::InputTag(label, CalibPrefix_ + "rawid" + CalibSuffix_));
0433   localdirx_token_ = consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "localdirx" + CalibSuffix_));
0434   localdiry_token_ = consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "localdiry" + CalibSuffix_));
0435   localdirz_token_ = consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "localdirz" + CalibSuffix_));
0436   firststrip_token_ =
0437       consumes<std::vector<unsigned short>>(edm::InputTag(label, CalibPrefix_ + "firststrip" + CalibSuffix_));
0438   nstrips_token_ = consumes<std::vector<unsigned short>>(edm::InputTag(label, CalibPrefix_ + "nstrips" + CalibSuffix_));
0439   saturation_token_ = consumes<std::vector<bool>>(edm::InputTag(label, CalibPrefix_ + "saturation" + CalibSuffix_));
0440   overlapping_token_ = consumes<std::vector<bool>>(edm::InputTag(label, CalibPrefix_ + "overlapping" + CalibSuffix_));
0441   farfromedge_token_ = consumes<std::vector<bool>>(edm::InputTag(label, CalibPrefix_ + "farfromedge" + CalibSuffix_));
0442   charge_token_ = consumes<std::vector<unsigned int>>(edm::InputTag(label, CalibPrefix_ + "charge" + CalibSuffix_));
0443   path_token_ = consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "path" + CalibSuffix_));
0444   chargeoverpath_token_ =
0445       consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "chargeoverpath" + CalibSuffix_));
0446   amplitude_token_ =
0447       consumes<std::vector<unsigned char>>(edm::InputTag(label, CalibPrefix_ + "amplitude" + CalibSuffix_));
0448   gainused_token_ = consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "gainused" + CalibSuffix_));
0449   gainusedTick_token_ =
0450       consumes<std::vector<double>>(edm::InputTag(label, CalibPrefix_ + "gainusedTick" + CalibSuffix_));
0451 
0452   edm::ParameterSet evtinfo_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("evtinfo");
0453   label = evtinfo_pset.getUntrackedParameter<string>("label");
0454   EventPrefix_ = evtinfo_pset.getUntrackedParameter<string>("prefix");
0455   EventSuffix_ = evtinfo_pset.getUntrackedParameter<string>("suffix");
0456   TrigTech_token_ = consumes<std::vector<bool>>(edm::InputTag(label, EventPrefix_ + "TrigTech" + EventSuffix_));
0457 
0458   edm::ParameterSet track_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("tracks");
0459   label = track_pset.getUntrackedParameter<string>("label");
0460   TrackPrefix_ = track_pset.getUntrackedParameter<string>("prefix");
0461   TrackSuffix_ = track_pset.getUntrackedParameter<string>("suffix");
0462 
0463   trackchi2ndof_token_ = consumes<std::vector<double>>(edm::InputTag(label, TrackPrefix_ + "chi2ndof" + TrackSuffix_));
0464   trackp_token_ = consumes<std::vector<float>>(edm::InputTag(label, TrackPrefix_ + "momentum" + TrackSuffix_));
0465   trackpt_token_ = consumes<std::vector<float>>(edm::InputTag(label, TrackPrefix_ + "pt" + TrackSuffix_));
0466   tracketa_token_ = consumes<std::vector<double>>(edm::InputTag(label, TrackPrefix_ + "eta" + TrackSuffix_));
0467   trackphi_token_ = consumes<std::vector<double>>(edm::InputTag(label, TrackPrefix_ + "phi" + TrackSuffix_));
0468   trackhitsvalid_token_ =
0469       consumes<std::vector<unsigned int>>(edm::InputTag(label, TrackPrefix_ + "hitsvalid" + TrackSuffix_));
0470   trackalgo_token_ = consumes<std::vector<int>>(edm::InputTag(label, TrackPrefix_ + "algo" + TrackSuffix_));
0471 
0472   tTopoToken_ = esConsumes<edm::Transition::BeginRun>();
0473   tkGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0474   runInfoToken_ = esConsumes<edm::Transition::BeginRun>();
0475   gainToken_ = esConsumes<edm::Transition::BeginRun>();
0476   qualityToken_ = esConsumes<edm::Transition::BeginRun>();
0477 }
0478 
0479 void SiStripGainFromCalibTree::bookDQMHistos(const char* dqm_dir, const char* tag) {
0480   edm::LogInfo("SiStripGainFromCalibTree")
0481       << "Setting " << dqm_dir << "in DQM and booking histograms for tag " << tag << std::endl;
0482 
0483   if (strcmp(booked_dir_.c_str(), dqm_dir) != 0) {
0484     booked_dir_ = dqm_dir;
0485     dbe->setCurrentFolder(dqm_dir);
0486   }
0487 
0488   std::string stag(tag);
0489   if (!stag.empty() && stag[0] != '_')
0490     stag.insert(0, 1, '_');
0491 
0492   std::string cvi = std::string("Charge_Vs_Index") + stag;
0493   //std::string cviA     = std::string("Charge_Vs_Index_Absolute")  + stag;
0494   std::string cvpTIB = std::string("Charge_Vs_PathlengthTIB") + stag;
0495   std::string cvpTOB = std::string("Charge_Vs_PathlengthTOB") + stag;
0496   std::string cvpTIDP = std::string("Charge_Vs_PathlengthTIDP") + stag;
0497   std::string cvpTIDM = std::string("Charge_Vs_PathlengthTIDM") + stag;
0498   std::string cvpTECP1 = std::string("Charge_Vs_PathlengthTECP1") + stag;
0499   std::string cvpTECP2 = std::string("Charge_Vs_PathlengthTECP2") + stag;
0500   std::string cvpTECM1 = std::string("Charge_Vs_PathlengthTECM1") + stag;
0501   std::string cvpTECM2 = std::string("Charge_Vs_PathlengthTECM2") + stag;
0502 
0503   int elepos = (m_harvestingMode && AlgoMode == "PCL") ? Harvest : statCollectionFromMode(tag);
0504 
0505   // The cluster charge is stored by exploiting a non uniform binning in order
0506   // reduce the histogram memory size. The bin width is relaxed with a falling
0507   // exponential function and the bin boundaries are stored in the binYarray.
0508   // The binXarray is used to provide as many bins as the APVs.
0509   //
0510   // More details about this implementations are here:
0511   // https://indico.cern.ch/event/649344/contributions/2672267/attachments/1498323/2332518/OptimizeChHisto.pdf
0512 
0513   std::vector<float> binXarray;
0514   binXarray.reserve(NStripAPVs + 1);
0515   for (int a = 0; a <= NStripAPVs; a++) {
0516     binXarray.push_back((float)a);
0517   }
0518 
0519   std::array<float, 688> binYarray;
0520   double p0 = 5.445;
0521   double p1 = 0.002113;
0522   double p2 = 69.01576;
0523   double y = 0.;
0524   for (int b = 0; b < 687; b++) {
0525     binYarray[b] = y;
0526     if (y <= 902.)
0527       y = y + 2.;
0528     else
0529       y = (p0 - log(exp(p0 - p1 * y) - p2 * p1)) / p1;
0530   }
0531   binYarray[687] = 4000.;
0532 
0533   Charge_Vs_Index[elepos] = dbe->book2S(cvi.c_str(), cvi.c_str(), NStripAPVs, &binXarray[0], 687, binYarray.data());
0534   //Charge_Vs_Index_Absolute[elepos]  = dbe->book2S(cviA.c_str()    , cviA.c_str()    , 88625, 0   , 88624,1000,0,4000);
0535   Charge_Vs_PathlengthTIB[elepos] = dbe->book2S(cvpTIB.c_str(), cvpTIB.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0536   Charge_Vs_PathlengthTOB[elepos] = dbe->book2S(cvpTOB.c_str(), cvpTOB.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0537   Charge_Vs_PathlengthTIDP[elepos] = dbe->book2S(cvpTIDP.c_str(), cvpTIDP.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0538   Charge_Vs_PathlengthTIDM[elepos] = dbe->book2S(cvpTIDM.c_str(), cvpTIDM.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0539   Charge_Vs_PathlengthTECP1[elepos] = dbe->book2S(cvpTECP1.c_str(), cvpTECP1.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0540   Charge_Vs_PathlengthTECP2[elepos] = dbe->book2S(cvpTECP2.c_str(), cvpTECP2.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0541   Charge_Vs_PathlengthTECM1[elepos] = dbe->book2S(cvpTECM1.c_str(), cvpTECM1.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0542   Charge_Vs_PathlengthTECM2[elepos] = dbe->book2S(cvpTECM2.c_str(), cvpTECM2.c_str(), 20, 0.3, 1.3, 250, 0, 2000);
0543 
0544   //Book Charge monitoring histograms
0545   std::vector<std::pair<std::string, std::string>> hnames =
0546       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0547   for (unsigned int i = 0; i < hnames.size(); i++) {
0548     std::string htag = (hnames[i]).first + stag;
0549     MonitorElement* monitor = dbe->book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.);
0550     int thick = APVGain::thickness((hnames[i]).first);
0551     int id = APVGain::subdetectorId((hnames[i]).first);
0552     int side = APVGain::subdetectorSide((hnames[i]).first);
0553     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0554     Charge_1[elepos].push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0555   }
0556 
0557   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG2");
0558   for (unsigned int i = 0; i < hnames.size(); i++) {
0559     std::string htag = (hnames[i]).first + stag;
0560     MonitorElement* monitor = dbe->book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.);
0561     int thick = APVGain::thickness((hnames[i]).first);
0562     int id = APVGain::subdetectorId((hnames[i]).first);
0563     int side = APVGain::subdetectorSide((hnames[i]).first);
0564     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0565     Charge_2[elepos].push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0566   }
0567 
0568   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1");
0569   for (unsigned int i = 0; i < hnames.size(); i++) {
0570     std::string htag = (hnames[i]).first + stag;
0571     MonitorElement* monitor = dbe->book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.);
0572     int thick = APVGain::thickness((hnames[i]).first);
0573     int id = APVGain::subdetectorId((hnames[i]).first);
0574     int side = APVGain::subdetectorSide((hnames[i]).first);
0575     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0576     Charge_3[elepos].push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0577   }
0578 
0579   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1G2");
0580   for (unsigned int i = 0; i < hnames.size(); i++) {
0581     std::string htag = (hnames[i]).first + stag;
0582     MonitorElement* monitor = dbe->book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.);
0583     int thick = APVGain::thickness((hnames[i]).first);
0584     int id = APVGain::subdetectorId((hnames[i]).first);
0585     int side = APVGain::subdetectorSide((hnames[i]).first);
0586     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0587     Charge_4[elepos].push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0588   }
0589 
0590   //Book validation histograms
0591   if (m_harvestingMode) {
0592     int MPVbin = 300;
0593     float MPVmin = 0.;
0594     float MPVmax = 600.;
0595 
0596     MPV_Vs_EtaTIB = dbe->book2DD("MPV_vs_EtaTIB", "MPV vs Eta TIB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0597     MPV_Vs_EtaTID = dbe->book2DD("MPV_vs_EtaTID", "MPV vs Eta TID", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0598     MPV_Vs_EtaTOB = dbe->book2DD("MPV_vs_EtaTOB", "MPV vs Eta TOB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0599     MPV_Vs_EtaTEC = dbe->book2DD("MPV_vs_EtaTEC", "MPV vs Eta TEC", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0600     MPV_Vs_EtaTECthin = dbe->book2DD("MPV_vs_EtaTEC1", "MPV vs Eta TEC-thin", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0601     MPV_Vs_EtaTECthick = dbe->book2DD("MPV_vs_EtaTEC2", "MPV vs Eta TEC-thick", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0602 
0603     MPV_Vs_PhiTIB = dbe->book2DD("MPV_vs_PhiTIB", "MPV vs Phi TIB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0604     MPV_Vs_PhiTID = dbe->book2DD("MPV_vs_PhiTID", "MPV vs Phi TID", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0605     MPV_Vs_PhiTOB = dbe->book2DD("MPV_vs_PhiTOB", "MPV vs Phi TOB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0606     MPV_Vs_PhiTEC = dbe->book2DD("MPV_vs_PhiTEC", "MPV vs Phi TEC", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0607     MPV_Vs_PhiTECthin = dbe->book2DD("MPV_vs_PhiTEC1", "MPV vs Phi TEC-thin", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0608     MPV_Vs_PhiTECthick = dbe->book2DD("MPV_vs_PhiTEC2", "MPV vs Phi TEC-thick", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0609 
0610     NoMPVfit = dbe->book2DD("NoMPVfit", "Modules with bad Landau Fit", 350, -350, 350, 240, 0, 120);
0611     NoMPVmasked = dbe->book2DD("NoMPVmasked", "Masked Modules", 350, -350, 350, 240, 0, 120);
0612 
0613     Gains = dbe->book1DD("Gains", "Gains", 300, 0, 2);
0614     MPVs = dbe->book1DD("MPVs", "MPVs", MPVbin, MPVmin, MPVmax);
0615     MPVs320 = dbe->book1DD("MPV_320", "MPV 320 thickness", MPVbin, MPVmin, MPVmax);
0616     MPVs500 = dbe->book1DD("MPV_500", "MPV 500 thickness", MPVbin, MPVmin, MPVmax);
0617     MPVsTIB = dbe->book1DD("MPV_TIB", "MPV TIB", MPVbin, MPVmin, MPVmax);
0618     MPVsTID = dbe->book1DD("MPV_TID", "MPV TID", MPVbin, MPVmin, MPVmax);
0619     MPVsTIDP = dbe->book1DD("MPV_TIDP", "MPV TIDP", MPVbin, MPVmin, MPVmax);
0620     MPVsTIDM = dbe->book1DD("MPV_TIDM", "MPV TIDM", MPVbin, MPVmin, MPVmax);
0621     MPVsTOB = dbe->book1DD("MPV_TOB", "MPV TOB", MPVbin, MPVmin, MPVmax);
0622     MPVsTEC = dbe->book1DD("MPV_TEC", "MPV TEC", MPVbin, MPVmin, MPVmax);
0623     MPVsTECP = dbe->book1DD("MPV_TECP", "MPV TECP", MPVbin, MPVmin, MPVmax);
0624     MPVsTECM = dbe->book1DD("MPV_TECM", "MPV TECM", MPVbin, MPVmin, MPVmax);
0625     MPVsTECthin = dbe->book1DD("MPV_TEC1", "MPV TEC1", MPVbin, MPVmin, MPVmax);
0626     MPVsTECthick = dbe->book1DD("MPV_TEC2", "MPV TEC2", MPVbin, MPVmin, MPVmax);
0627     MPVsTECP1 = dbe->book1DD("MPV_TECP1", "MPV TECP1", MPVbin, MPVmin, MPVmax);
0628     MPVsTECP2 = dbe->book1DD("MPV_TECP2", "MPV TECP2", MPVbin, MPVmin, MPVmax);
0629     MPVsTECM1 = dbe->book1DD("MPV_TECM1", "MPV TECM1", MPVbin, MPVmin, MPVmax);
0630     MPVsTECM2 = dbe->book1DD("MPV_TECM2", "MPV TECM2", MPVbin, MPVmin, MPVmax);
0631 
0632     MPVError = dbe->book1DD("MPVError", "MPV Error", 150, 0, 150);
0633     MPVErrorVsMPV = dbe->book2DD("MPVErrorVsMPV", "MPV Error vs MPV", 300, 0, 600, 150, 0, 150);
0634     MPVErrorVsEta = dbe->book2DD("MPVErrorVsEta", "MPV Error vs Eta", 50, -3.0, 3.0, 150, 0, 150);
0635     MPVErrorVsPhi = dbe->book2DD("MPVErrorVsPhi", "MPV Error vs Phi", 50, -3.4, 3.4, 150, 0, 150);
0636     MPVErrorVsN = dbe->book2DD("MPVErrorVsN", "MPV Error vs N", 500, 0, 1000, 150, 0, 150);
0637 
0638     DiffWRTPrevGainTIB = dbe->book1DD("DiffWRTPrevGainTIB", "Diff w.r.t. PrevGain TIB", 250, 0, 2);
0639     DiffWRTPrevGainTID = dbe->book1DD("DiffWRTPrevGainTID", "Diff w.r.t. PrevGain TID", 250, 0, 2);
0640     DiffWRTPrevGainTOB = dbe->book1DD("DiffWRTPrevGainTOB", "Diff w.r.t. PrevGain TOB", 250, 0, 2);
0641     DiffWRTPrevGainTEC = dbe->book1DD("DiffWRTPrevGainTEC", "Diff w.r.t. PrevGain TEC", 250, 0, 2);
0642 
0643     GainVsPrevGainTIB = dbe->book2DD("GainVsPrevGainTIB", "Gain vs PrevGain TIB", 100, 0, 2, 100, 0, 2);
0644     GainVsPrevGainTID = dbe->book2DD("GainVsPrevGainTID", "Gain vs PrevGain TID", 100, 0, 2, 100, 0, 2);
0645     GainVsPrevGainTOB = dbe->book2DD("GainVsPrevGainTOB", "Gain vs PrevGain TOB", 100, 0, 2, 100, 0, 2);
0646     GainVsPrevGainTEC = dbe->book2DD("GainVsPrevGainTEC", "Gain vs PrevGain TEC", 100, 0, 2, 100, 0, 2);
0647 
0648     std::vector<std::pair<std::string, std::string>> hnames =
0649         APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "newG2");
0650     for (unsigned int i = 0; i < hnames.size(); i++) {
0651       MonitorElement* monitor = dbe->book1DD((hnames[i]).first.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.);
0652       int thick = APVGain::thickness((hnames[i]).first);
0653       int id = APVGain::subdetectorId((hnames[i]).first);
0654       int side = APVGain::subdetectorSide((hnames[i]).first);
0655       int plane = APVGain::subdetectorPlane((hnames[i]).first);
0656       newCharge.push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0657     }
0658   }
0659 }
0660 
0661 void SiStripGainFromCalibTree::algoBeginJob(const edm::EventSetup& iSetup) {
0662   edm::LogInfo("SiStripGainFromCalibTree") << "AlgoMode        : " << AlgoMode << "\n"
0663                                            << "CalibrationMode : " << m_calibrationMode << "\n"
0664                                            << "HarvestingMode  : " << m_harvestingMode << std::endl;
0665   //Setup DQM histograms
0666   if (AlgoMode != "PCL" or m_harvestingMode) {
0667     const char* dqm_dir = "AlCaReco/SiStripGainsHarvesting/";
0668     this->bookDQMHistos(dqm_dir, dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str());
0669   } else {
0670     //Check consistency of calibration Mode and BField only for the ALCAPROMPT in the PCL workflow
0671     if (!isBFieldConsistentWithMode(iSetup)) {
0672       string prevMode = m_calibrationMode;
0673       swapBFieldMode();
0674       edm::LogInfo("SiStripGainFromCalibTree") << "Switching calibration mode for endorsing BField status: " << prevMode
0675                                                << " ==> " << m_calibrationMode << std::endl;
0676     }
0677     std::string dqm_dir = m_DQMdir + ((m_splitDQMstat) ? m_calibrationMode : "") + "/";
0678     int elem = statCollectionFromMode(m_calibrationMode.c_str());
0679     this->bookDQMHistos(dqm_dir.c_str(), dqm_tag_[elem].c_str());
0680   }
0681 
0682   tTopo_ = &iSetup.getData(tTopoToken_);
0683 
0684   auto const& Det = iSetup.getData(tkGeomToken_).dets();
0685 
0686   NPixelDets = 0;
0687   NStripAPVs = 0;
0688   unsigned int Index = 0;
0689   for (unsigned int i = 0; i < Det.size(); i++) {
0690     DetId Detid = Det[i]->geographicalId();
0691     int SubDet = Detid.subdetId();
0692 
0693     if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
0694         SubDet == StripSubdetector::TEC) {
0695       auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
0696       if (!DetUnit)
0697         continue;
0698 
0699       const StripTopology& Topo = DetUnit->specificTopology();
0700       unsigned int NAPV = Topo.nstrips() / 128;
0701 
0702       for (unsigned int j = 0; j < NAPV; j++) {
0703         stAPVGain* APV = new stAPVGain;
0704         APV->Index = Index;
0705         APV->Bin = -1;
0706         APV->DetId = Detid.rawId();
0707         APV->APVId = j;
0708         APV->SubDet = SubDet;
0709         APV->FitMPV = -1;
0710         APV->FitMPVErr = -1;
0711         APV->FitWidth = -1;
0712         APV->FitWidthErr = -1;
0713         APV->FitChi2 = -1;
0714         APV->FitNorm = -1;
0715         APV->Gain = -1;
0716         APV->PreviousGain = 1;
0717         APV->PreviousGainTick = 1;
0718         APV->x = DetUnit->position().basicVector().x();
0719         APV->y = DetUnit->position().basicVector().y();
0720         APV->z = DetUnit->position().basicVector().z();
0721         APV->Eta = DetUnit->position().basicVector().eta();
0722         APV->Phi = DetUnit->position().basicVector().phi();
0723         APV->R = DetUnit->position().basicVector().transverse();
0724         APV->Thickness = DetUnit->surface().bounds().thickness();
0725         APV->NEntries = 0;
0726         APV->isMasked = false;
0727 
0728         APVsCollOrdered.push_back(APV);
0729         APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0730         Index++;
0731         NStripAPVs++;
0732       }
0733     }
0734   }
0735 
0736   for (unsigned int i = 0; i < Det.size();
0737        i++) {  //Make two loop such that the Pixel information is added at the end --> make transition simpler
0738     DetId Detid = Det[i]->geographicalId();
0739     int SubDet = Detid.subdetId();
0740     if (SubDet == PixelSubdetector::PixelBarrel || SubDet == PixelSubdetector::PixelEndcap) {
0741       auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
0742       if (!DetUnit)
0743         continue;
0744 
0745       const PixelTopology& Topo = DetUnit->specificTopology();
0746       unsigned int NROCRow = Topo.nrows() / (80.);
0747       unsigned int NROCCol = Topo.ncolumns() / (52.);
0748 
0749       for (unsigned int j = 0; j < NROCRow; j++) {
0750         for (unsigned int i = 0; i < NROCCol; i++) {
0751           stAPVGain* APV = new stAPVGain;
0752           APV->Index = Index;
0753           APV->Bin = -1;
0754           APV->DetId = Detid.rawId();
0755           APV->APVId = (j << 3 | i);
0756           APV->SubDet = SubDet;
0757           APV->FitMPV = -1;
0758           APV->FitMPVErr = -1;
0759           APV->FitWidth = -1;
0760           APV->FitWidthErr = -1;
0761           APV->FitChi2 = -1;
0762           APV->Gain = -1;
0763           APV->PreviousGain = 1;
0764           APV->PreviousGainTick = 1;
0765           APV->x = DetUnit->position().basicVector().x();
0766           APV->y = DetUnit->position().basicVector().y();
0767           APV->z = DetUnit->position().basicVector().z();
0768           APV->Eta = DetUnit->position().basicVector().eta();
0769           APV->Phi = DetUnit->position().basicVector().phi();
0770           APV->R = DetUnit->position().basicVector().transverse();
0771           APV->Thickness = DetUnit->surface().bounds().thickness();
0772           APV->isMasked = false;  //SiPixelQuality_->IsModuleBad(Detid.rawId());
0773           APV->NEntries = 0;
0774 
0775           APVsCollOrdered.push_back(APV);
0776           APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0777           Index++;
0778           NPixelDets++;
0779         }
0780       }
0781     }
0782   }
0783 
0784   MakeCalibrationMap();
0785 
0786   NEvent = 0;
0787   NTrack = 0;
0788   NClusterStrip = 0;
0789   NClusterPixel = 0;
0790   SRun = 1 << 31;
0791   ERun = 0;
0792   GOOD = 0;
0793   BAD = 0;
0794   MASKED = 0;
0795 }
0796 
0797 bool SiStripGainFromCalibTree::isBFieldConsistentWithMode(const edm::EventSetup& iSetup) const {
0798   const auto& runInfo = iSetup.getData(runInfoToken_);
0799 
0800   double average_current = runInfo.m_avg_current;
0801   bool isOn = (average_current > MagFieldCurrentTh);
0802   bool is0T = (m_calibrationMode.substr(m_calibrationMode.length() - 2) == "0T");
0803 
0804   return ((isOn && !is0T) || (!isOn && is0T));
0805 }
0806 
0807 void SiStripGainFromCalibTree::swapBFieldMode() {
0808   if (m_calibrationMode.substr(m_calibrationMode.length() - 2) == "0T") {
0809     m_calibrationMode.erase(m_calibrationMode.length() - 2, 2);
0810   } else {
0811     m_calibrationMode.append("0T");
0812   }
0813 }
0814 
0815 void SiStripGainFromCalibTree::algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0816   if (!m_harvestingMode && AlgoMode == "PCL") {
0817     //Check consistency of calibration Mode and BField only for the ALCAPROMPT in the PCL workflow
0818     if (!isBFieldConsistentWithMode(iSetup)) {
0819       string prevMode = m_calibrationMode;
0820       swapBFieldMode();
0821       edm::LogInfo("SiStripGainFromCalibTree") << "Switching calibration mode for endorsing BField status: " << prevMode
0822                                                << " ==> " << m_calibrationMode << std::endl;
0823     }
0824   }
0825 
0826   const auto gainHandle = iSetup.getHandle(gainToken_);
0827   if (!gainHandle.isValid()) {
0828     edm::LogError("SiStripGainFromCalibTree") << "gainHandle is not valid\n";
0829     exit(0);
0830   }
0831 
0832   const auto& siStripQuality = iSetup.getData(qualityToken_);
0833   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0834     stAPVGain* APV = APVsCollOrdered[a];
0835 
0836     // MM. 03/03/2017 all of this makes sense only for SiStrip (i.e. get me out of here if I am on a pixel ROC)
0837     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0838       continue;
0839 
0840     APV->isMasked = siStripQuality.IsApvBad(APV->DetId, APV->APVId);
0841     //      if(!FirstSetOfConstants){
0842 
0843     if (gainHandle->getNumberOfTags() != 2) {
0844       edm::LogError("SiStripGainFromCalibTree") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
0845       fflush(stdout);
0846       exit(0);
0847     };
0848     float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
0849     if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
0850       edm::LogWarning("SiStripGainFromCalibTree") << "WARNING: ParticleGain in the global tag changed\n";
0851     APV->PreviousGain = newPreviousGain;
0852 
0853     float newPreviousGainTick = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
0854     if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
0855       edm::LogWarning("SiStripGainFromCalibTree")
0856           << "WARNING: TickMarkGain in the global tag changed\n"
0857           << std::endl
0858           << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
0859           << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
0860           << std::endl;
0861     }
0862     APV->PreviousGainTick = newPreviousGainTick;
0863 
0864     //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f (G1) x %8.4f (G2)\n",APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain);
0865     //      }
0866   }
0867 }
0868 
0869 void SiStripGainFromCalibTree::algoEndRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0870   if (AlgoMode == "PCL" && !m_harvestingMode)
0871     return;  //nothing to do in that case
0872 
0873   if (AlgoMode == "PCL" and m_harvestingMode) {
0874     // Load the 2D histograms from the DQM objects
0875     // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
0876     // DQMStore so that they can be used in the fit
0877 
0878     edm::LogInfo("SiStripGainFromCalibTree") << "Starting harvesting statistics" << std::endl;
0879 
0880     // check the required tag before adding histograms
0881     int elepos = statCollectionFromMode(m_calibrationMode.c_str());
0882     if (elepos != Harvest) {
0883       //collect statistics from DQM into the related monitored elements
0884       std::string stag = m_calibrationMode;
0885       if (!stag.empty() && stag[0] != '_')
0886         stag.insert(0, 1, '_');
0887 
0888       if (elepos == -1) {
0889         //implememt backward compatibility
0890         elepos = 0;
0891         stag = "";
0892       }
0893 
0894       std::string DQM_dir = m_DQMdir + ((m_splitDQMstat) ? m_calibrationMode : "");
0895 
0896       std::string cvi = DQM_dir + std::string("/Charge_Vs_Index") + stag;
0897       //std::string cviA     = DQM_dir + std::string("/Charge_Vs_Index_Absolute")  + stag;
0898       std::string cvpTIB = DQM_dir + std::string("/Charge_Vs_PathlengthTIB") + stag;
0899       std::string cvpTOB = DQM_dir + std::string("/Charge_Vs_PathlengthTOB") + stag;
0900       std::string cvpTIDP = DQM_dir + std::string("/Charge_Vs_PathlengthTIDP") + stag;
0901       std::string cvpTIDM = DQM_dir + std::string("/Charge_Vs_PathlengthTIDM") + stag;
0902       std::string cvpTECP1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP1") + stag;
0903       std::string cvpTECP2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP2") + stag;
0904       std::string cvpTECM1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM1") + stag;
0905       std::string cvpTECM2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM2") + stag;
0906 
0907       Charge_Vs_Index[elepos] = dbe->get(cvi);
0908       //Charge_Vs_Index_Absolute[elepos]  = dbe->get(cviA.c_str());
0909       Charge_Vs_PathlengthTIB[elepos] = dbe->get(cvpTIB);
0910       Charge_Vs_PathlengthTOB[elepos] = dbe->get(cvpTOB);
0911       Charge_Vs_PathlengthTIDP[elepos] = dbe->get(cvpTIDP);
0912       Charge_Vs_PathlengthTIDM[elepos] = dbe->get(cvpTIDM);
0913       Charge_Vs_PathlengthTECP1[elepos] = dbe->get(cvpTECP1);
0914       Charge_Vs_PathlengthTECP2[elepos] = dbe->get(cvpTECP2);
0915       Charge_Vs_PathlengthTECM1[elepos] = dbe->get(cvpTECM1);
0916       Charge_Vs_PathlengthTECM2[elepos] = dbe->get(cvpTECM2);
0917 
0918       if (Charge_Vs_Index[elepos] == nullptr) {
0919         edm::LogError("SiStripGainFromCalibTree")
0920             << "Harvesting: could not retrieve " << cvi.c_str() << ", statistics will not be summed!" << std::endl;
0921       } else {
0922         merge((Charge_Vs_Index[Harvest])->getTH2S(), (Charge_Vs_Index[elepos])->getTH2S());
0923         edm::LogInfo("SiStripGainFromCalibTree")
0924             << "Harvesting " << (Charge_Vs_Index[elepos])->getTH2S()->GetEntries() << " more clusters" << std::endl;
0925       }
0926 
0927       //if (Charge_Vs_Index_Absolute[elepos]==0) {
0928       //    edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cviA.c_str()
0929       //                                              << ", statistics will not be summed!" << std::endl;
0930       //} else merge( (Charge_Vs_Index_Absolute[Harvest])->getTH2S(), (Charge_Vs_Index_Absolute[elepos])->getTH2S() );
0931 
0932       if (Charge_Vs_PathlengthTIB[elepos] == nullptr) {
0933         edm::LogError("SiStripGainFromCalibTree")
0934             << "Harvesting: could not retrieve " << cvpTIB.c_str() << ", statistics will not be summed!" << std::endl;
0935       } else
0936         (Charge_Vs_PathlengthTIB[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIB[elepos])->getTH2S());
0937 
0938       if (Charge_Vs_PathlengthTOB[elepos] == nullptr) {
0939         edm::LogError("SiStripGainFromCalibTree")
0940             << "Harvesting: could not retrieve " << cvpTOB.c_str() << ", statistics will not be summed!" << std::endl;
0941       } else
0942         (Charge_Vs_PathlengthTOB[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTOB[elepos])->getTH2S());
0943 
0944       if (Charge_Vs_PathlengthTIDP[elepos] == nullptr) {
0945         edm::LogError("SiStripGainFromCalibTree")
0946             << "Harvesting: could not retrieve " << cvpTIDP.c_str() << ", statistics will not be summed!" << std::endl;
0947       } else
0948         (Charge_Vs_PathlengthTIDP[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIDP[elepos])->getTH2S());
0949 
0950       if (Charge_Vs_PathlengthTIDM[elepos] == nullptr) {
0951         edm::LogError("SiStripGainFromCalibTree")
0952             << "Harvesting: could not retrieve " << cvpTIDM.c_str() << ", statistics will not be summed!" << std::endl;
0953       } else
0954         (Charge_Vs_PathlengthTIDM[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIDM[elepos])->getTH2S());
0955 
0956       if (Charge_Vs_PathlengthTECP1[elepos] == nullptr) {
0957         edm::LogError("SiStripGainFromCalibTree")
0958             << "Harvesting: could not retrieve " << cvpTECP1.c_str() << ", statistics will not be summed!" << std::endl;
0959       } else
0960         (Charge_Vs_PathlengthTECP1[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECP1[elepos])->getTH2S());
0961 
0962       if (Charge_Vs_PathlengthTECP2[elepos] == nullptr) {
0963         edm::LogError("SiStripGainFromCalibTree")
0964             << "Harvesting: could not retrieve " << cvpTECP2.c_str() << ", statistics will not be summed!" << std::endl;
0965       } else
0966         (Charge_Vs_PathlengthTECP2[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECP2[elepos])->getTH2S());
0967 
0968       if (Charge_Vs_PathlengthTECM1[elepos] == nullptr) {
0969         edm::LogError("SiStripGainFromCalibTree")
0970             << "Harvesting: could not retrieve " << cvpTECM1.c_str() << ", statistics will not be summed!" << std::endl;
0971       } else
0972         (Charge_Vs_PathlengthTECM1[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECM1[elepos])->getTH2S());
0973 
0974       if (Charge_Vs_PathlengthTECM2[elepos] == nullptr) {
0975         edm::LogError("SiStripGainFromCalibTree")
0976             << "Harvesting: could not retrieve " << cvpTECM2.c_str() << ", statistics will not be summed!" << std::endl;
0977       } else
0978         (Charge_Vs_PathlengthTECM2[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECM2[elepos])->getTH2S());
0979 
0980       // Gather Charge monitoring histograms
0981       std::vector<std::pair<std::string, std::string>> tags =
0982           APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0983       for (unsigned int i = 0; i < tags.size(); i++) {
0984         std::string tag = DQM_dir + "/" + (tags[i]).first + stag;
0985         Charge_1[elepos].push_back(APVGain::APVmon(0, 0, 0, 0, dbe->get(tag)));
0986         if ((Charge_1[elepos][i]).getMonitor() == nullptr) {
0987           edm::LogError("SiStripGainFromCalibTree")
0988               << "Harvesting: could not retrieve " << tag.c_str() << ", statistics will not be summed!" << std::endl;
0989         } else
0990           (Charge_1[Harvest][i]).getMonitor()->getTH1D()->Add((Charge_1[elepos][i]).getMonitor()->getTH1D());
0991       }
0992 
0993       tags = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG2");
0994       for (unsigned int i = 0; i < tags.size(); i++) {
0995         std::string tag = DQM_dir + "/" + (tags[i]).first + stag;
0996         Charge_2[elepos].push_back(APVGain::APVmon(0, 0, 0, 0, dbe->get(tag)));
0997         if ((Charge_2[elepos][i]).getMonitor() == nullptr) {
0998           edm::LogError("SiStripGainFromCalibTree")
0999               << "Harvesting: could not retrieve " << tag.c_str() << ", statistics will not be summed!" << std::endl;
1000         } else
1001           (Charge_2[Harvest][i]).getMonitor()->getTH1D()->Add((Charge_2[elepos][i]).getMonitor()->getTH1D());
1002       }
1003 
1004       tags = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1");
1005       for (unsigned int i = 0; i < tags.size(); i++) {
1006         std::string tag = DQM_dir + "/" + (tags[i]).first + stag;
1007         Charge_3[elepos].push_back(APVGain::APVmon(0, 0, 0, 0, dbe->get(tag)));
1008         if ((Charge_3[elepos][i]).getMonitor() == nullptr) {
1009           edm::LogError("SiStripGainFromCalibTree")
1010               << "Harvesting: could not retrieve " << tag.c_str() << ", statistics will not be summed!" << std::endl;
1011         } else
1012           (Charge_3[Harvest][i]).getMonitor()->getTH1D()->Add((Charge_3[elepos][i]).getMonitor()->getTH1D());
1013       }
1014 
1015       tags = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1G2");
1016       for (unsigned int i = 0; i < tags.size(); i++) {
1017         std::string tag = DQM_dir + "/" + (tags[i]).first + stag;
1018         Charge_4[elepos].push_back(APVGain::APVmon(0, 0, 0, 0, dbe->get(tag)));
1019         if ((Charge_4[elepos][i]).getMonitor() == nullptr) {
1020           edm::LogError("SiStripGainFromCalibTree")
1021               << "Harvesting: could not retrieve " << tag.c_str() << ", statistics will not be summed!" << std::endl;
1022         } else
1023           (Charge_4[Harvest][i]).getMonitor()->getTH1D()->Add((Charge_4[elepos][i]).getMonitor()->getTH1D());
1024       }
1025     }
1026   }
1027 }
1028 
1029 void SiStripGainFromCalibTree::algoEndJob() {
1030   if (AlgoMode == "PCL" && !m_harvestingMode)
1031     return;  //nothing to do in that case
1032 
1033   if (AlgoMode == "CalibTree") {
1034     edm::LogInfo("SiStripGainFromCalibTree") << "Analyzing calibration tree" << std::endl;
1035     // Loop on calibTrees to fill the 2D histograms
1036     algoAnalyzeTheTree();
1037   } else if (m_harvestingMode) {
1038     NClusterStrip = (Charge_Vs_Index[Harvest])->getTH2S()->Integral(0, NStripAPVs + 1, 0, 99999);
1039     //NClusterPixel = (Charge_Vs_Index[Harvest])->getTH2S()->Integral(NStripAPVs+2, NStripAPVs+NPixelDets+2, 0, 99999 );
1040   }
1041 
1042   // Now that we have the full statistics we can extract the information of the 2D histograms
1043   algoComputeMPVandGain();
1044 
1045   // Result monitoring
1046   qualityMonitor();
1047 
1048   // Force the DB object writing,
1049   // thus setting the IOV as the first processed run (if timeFromEndRun is set to false)
1050   storeOnDbNow();
1051 
1052   if (AlgoMode != "PCL" or saveSummary) {
1053     edm::LogInfo("SiStripGainFromCalibTree") << "Saving summary into root file" << std::endl;
1054 
1055     //also save the 2D monitor elements to this file as TH2D tfs
1056     tfs = edm::Service<TFileService>().operator->();
1057 
1058     //save only the statistics for the calibrationTag
1059     int elepos = statCollectionFromMode(m_calibrationMode.c_str());
1060 
1061     if (Charge_Vs_Index[elepos] != nullptr)
1062       tfs->make<TH2S>(*(Charge_Vs_Index[elepos])->getTH2S());
1063     //if( Charge_Vs_Index_Absolute[elepos]!=0 )  tfs->make<TH2S> ( *(Charge_Vs_Index_Absolute[elepos])->getTH2S() );
1064     if (Charge_Vs_PathlengthTIB[elepos] != nullptr)
1065       tfs->make<TH2S>(*(Charge_Vs_PathlengthTIB[elepos])->getTH2S());
1066     if (Charge_Vs_PathlengthTOB[elepos] != nullptr)
1067       tfs->make<TH2S>(*(Charge_Vs_PathlengthTOB[elepos])->getTH2S());
1068     if (Charge_Vs_PathlengthTIDP[elepos] != nullptr)
1069       tfs->make<TH2S>(*(Charge_Vs_PathlengthTIDP[elepos])->getTH2S());
1070     if (Charge_Vs_PathlengthTIDM[elepos] != nullptr)
1071       tfs->make<TH2S>(*(Charge_Vs_PathlengthTIDM[elepos])->getTH2S());
1072     if (Charge_Vs_PathlengthTECP1[elepos] != nullptr)
1073       tfs->make<TH2S>(*(Charge_Vs_PathlengthTECP1[elepos])->getTH2S());
1074     if (Charge_Vs_PathlengthTECP2[elepos] != nullptr)
1075       tfs->make<TH2S>(*(Charge_Vs_PathlengthTECP2[elepos])->getTH2S());
1076     if (Charge_Vs_PathlengthTECM1[elepos] != nullptr)
1077       tfs->make<TH2S>(*(Charge_Vs_PathlengthTECM1[elepos])->getTH2S());
1078     if (Charge_Vs_PathlengthTECM2[elepos] != nullptr)
1079       tfs->make<TH2S>(*(Charge_Vs_PathlengthTECM2[elepos])->getTH2S());
1080 
1081     storeOnTree(tfs);
1082   }
1083 }
1084 
1085 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange) {
1086   FitResults[0] = -0.5;  //MPV
1087   FitResults[1] = 0;     //MPV error
1088   FitResults[2] = -0.5;  //Width
1089   FitResults[3] = 0;     //Width error
1090   FitResults[4] = -0.5;  //Fit Chi2/NDF
1091   FitResults[5] = 0;     //Normalization
1092 
1093   if (InputHisto->GetEntries() < MinNrEntries)
1094     return;
1095 
1096   // perform fit with standard landau
1097   TF1* MyLandau = new TF1("MyLandau", "landau", LowRange, HighRange);
1098   MyLandau->SetParameter(1, 300);
1099   InputHisto->Fit(MyLandau, "0QR WW");
1100 
1101   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
1102   FitResults[0] = MyLandau->GetParameter(1);                      //MPV
1103   FitResults[1] = MyLandau->GetParError(1);                       //MPV error
1104   FitResults[2] = MyLandau->GetParameter(2);                      //Width
1105   FitResults[3] = MyLandau->GetParError(2);                       //Width error
1106   FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF();  //Fit Chi2/NDF
1107   FitResults[5] = MyLandau->GetParameter(0);
1108 
1109   delete MyLandau;
1110 }
1111 
1112 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults) {
1113   if (FitResults[0] <= 0)
1114     return false;
1115   //   if(FitResults[1] > MaxMPVError   )return false;
1116   //   if(FitResults[4] > MaxChi2OverNDF)return false;
1117   return true;
1118 }
1119 
1120 void SiStripGainFromCalibTree::processEvent() {
1121   edm::LogInfo("SiStripGainFromCalibTree") << "Processing run " << runnumber << " and event " << eventnumber << " for "
1122                                            << m_calibrationMode << " calibration." << std::endl;
1123 
1124   if (runnumber < SRun)
1125     SRun = runnumber;
1126   if (runnumber > ERun)
1127     ERun = runnumber;
1128 
1129   NEvent++;
1130   NTrack += (*trackp).size();
1131 
1132   int elepos = statCollectionFromMode(m_calibrationMode.c_str());
1133 
1134   unsigned int FirstAmplitude = 0;
1135   for (unsigned int i = 0; i < (*chargeoverpath).size(); i++) {
1136     FirstAmplitude += (*nstrips)[i];
1137     int TI = (*trackindex)[i];
1138 
1139     //printf("%i - %i - %i %i %i\n", (int)(*rawid)[i], (int)(*firststrip)[i]/128, (int)(*farfromedge)[i], (int)(*overlapping)[i], (int)(*saturation )[i] );
1140     if ((*tracketa)[TI] < MinTrackEta)
1141       continue;
1142     if ((*tracketa)[TI] > MaxTrackEta)
1143       continue;
1144     if ((*trackp)[TI] < MinTrackMomentum)
1145       continue;
1146     if ((*trackp)[TI] > MaxTrackMomentum)
1147       continue;
1148     if ((*trackhitsvalid)[TI] < MinTrackHits)
1149       continue;
1150     if ((*trackchi2ndof)[TI] > MaxTrackChiOverNdf)
1151       continue;
1152     if ((*trackalgo)[TI] > MaxTrackingIteration)
1153       continue;
1154 
1155     stAPVGain* APV =
1156         APVsColl[((*rawid)[i] << 4) |
1157                  ((*firststrip)[i] /
1158                   128)];  //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
1159 
1160     if (APV->SubDet > 2 && (*farfromedge)[i] == false)
1161       continue;
1162     if (APV->SubDet > 2 && (*overlapping)[i] == true)
1163       continue;
1164     if (APV->SubDet > 2 && (*saturation)[i] && !AllowSaturation)
1165       continue;
1166     if (APV->SubDet > 2 && (*nstrips)[i] > MaxNrStrips)
1167       continue;
1168 
1169     //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
1170 
1171     //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
1172     //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) +  pow((*localdirz)[i],2) ) ) * (180/3.14159265);
1173     //double beta  = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) +  pow((*localdirz)[i],2) ) ) * (180/3.14159265);
1174 
1175     //printf("NStrip = %i : Charge = %i --> Path = %f  --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
1176     //printf("Amplitudes: ");
1177     //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
1178     //printf("\n");
1179 
1180     if (APV->SubDet > 2) {
1181       NClusterStrip++;
1182     } else {
1183       NClusterPixel++;
1184     }
1185 
1186     int Charge = 0;
1187     if (APV->SubDet > 2 && (useCalibration || !FirstSetOfConstants)) {
1188       bool Saturation = false;
1189       for (unsigned int s = 0; s < (*nstrips)[i]; s++) {
1190         int StripCharge = (*amplitude)[FirstAmplitude - (*nstrips)[i] + s];
1191         if (useCalibration && !FirstSetOfConstants) {
1192           StripCharge = (int)(StripCharge * (APV->PreviousGain / APV->CalibGain));
1193         } else if (useCalibration) {
1194           StripCharge = (int)(StripCharge / APV->CalibGain);
1195         } else if (!FirstSetOfConstants) {
1196           StripCharge = (int)(StripCharge * APV->PreviousGain);
1197         }
1198         if (StripCharge > 1024) {
1199           StripCharge = 255;
1200           Saturation = true;
1201         } else if (StripCharge > 254) {
1202           StripCharge = 254;
1203           Saturation = true;
1204         }
1205         Charge += StripCharge;
1206       }
1207       if (Saturation && !AllowSaturation)
1208         continue;
1209     } else if (APV->SubDet > 2) {
1210       Charge = (*charge)[i];
1211     } else {
1212       Charge = (*charge)[i] / 265.0;  //expected scale factor between pixel and strip charge
1213     }
1214 
1215     //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
1216 
1217     double ClusterChargeOverPath = ((double)Charge) / (*path)[i];
1218     if (APV->SubDet > 2) {
1219       if (Validation) {
1220         ClusterChargeOverPath /= (*gainused)[i];
1221       }
1222       if (OldGainRemoving) {
1223         ClusterChargeOverPath *= (*gainused)[i];
1224       }
1225     }
1226 
1227     // keep pixel cluster charge processing until here
1228     if (APV->SubDet <= 2)
1229       continue;
1230 
1231     (Charge_Vs_Index[elepos])->Fill(APV->Index, ClusterChargeOverPath);
1232 
1233     // Compute the charge for monitoring and fill the relative histograms
1234     int mCharge1 = 0;
1235     int mCharge2 = 0;
1236     int mCharge3 = 0;
1237     int mCharge4 = 0;
1238     if (APV->SubDet > 2) {
1239       for (unsigned int s = 0; s < (*nstrips)[i]; s++) {
1240         int StripCharge = (*amplitude)[FirstAmplitude - (*nstrips)[i] + s];
1241         if (StripCharge > 1024)
1242           StripCharge = 255;
1243         else if (StripCharge > 254)
1244           StripCharge = 254;
1245         mCharge1 += StripCharge;
1246         mCharge2 += StripCharge;
1247         mCharge3 += StripCharge;
1248         mCharge4 += StripCharge;
1249       }
1250       // Revome gains for monitoring
1251       mCharge2 *= (*gainused)[i];                         // remove G2
1252       mCharge3 *= (*gainusedTick)[i];                     // remove G1
1253       mCharge4 *= ((*gainused)[i] * (*gainusedTick)[i]);  // remove G1 and G2
1254     }
1255     std::vector<APVGain::APVmon>& v1 = Charge_1[elepos];
1256     std::vector<MonitorElement*> cmon1 = APVGain::FetchMonitor(v1, (*rawid)[i], tTopo_);
1257     for (unsigned int m = 0; m < cmon1.size(); m++)
1258       cmon1[m]->Fill(((double)mCharge1) / (*path)[i]);
1259 
1260     std::vector<APVGain::APVmon>& v2 = Charge_2[elepos];
1261     std::vector<MonitorElement*> cmon2 = APVGain::FetchMonitor(v2, (*rawid)[i], tTopo_);
1262     for (unsigned int m = 0; m < cmon2.size(); m++)
1263       cmon2[m]->Fill(((double)mCharge2) / (*path)[i]);
1264 
1265     std::vector<APVGain::APVmon>& v3 = Charge_3[elepos];
1266     std::vector<MonitorElement*> cmon3 = APVGain::FetchMonitor(v3, (*rawid)[i], tTopo_);
1267     for (unsigned int m = 0; m < cmon3.size(); m++)
1268       cmon3[m]->Fill(((double)mCharge3) / (*path)[i]);
1269 
1270     std::vector<APVGain::APVmon>& v4 = Charge_4[elepos];
1271     std::vector<MonitorElement*> cmon4 = APVGain::FetchMonitor(v4, (*rawid)[i], tTopo_);
1272     for (unsigned int m = 0; m < cmon4.size(); m++)
1273       cmon4[m]->Fill(((double)mCharge4) / (*path)[i]);
1274 
1275     // Fill Charge Vs pathLenght histograms
1276     if (APV->SubDet == StripSubdetector::TIB) {
1277       (Charge_Vs_PathlengthTIB[elepos])->Fill((*path)[i], Charge);  // TIB
1278 
1279     } else if (APV->SubDet == StripSubdetector::TOB) {
1280       (Charge_Vs_PathlengthTOB[elepos])->Fill((*path)[i], Charge);  // TOB
1281 
1282     } else if (APV->SubDet == StripSubdetector::TID) {
1283       if (APV->Eta < 0) {
1284         (Charge_Vs_PathlengthTIDM[elepos])->Fill((*path)[i], Charge);
1285       }  // TID minus
1286       else if (APV->Eta > 0) {
1287         (Charge_Vs_PathlengthTIDP[elepos])->Fill((*path)[i], Charge);
1288       }  // TID plus
1289 
1290     } else if (APV->SubDet == StripSubdetector::TEC) {
1291       if (APV->Eta < 0) {
1292         if (APV->Thickness < 0.04) {
1293           (Charge_Vs_PathlengthTECM1[elepos])->Fill((*path)[i], Charge);
1294         }  // TEC minus, type 1
1295         else if (APV->Thickness > 0.04) {
1296           (Charge_Vs_PathlengthTECM2[elepos])->Fill((*path)[i], Charge);
1297         }  // TEC minus, type 2
1298       } else if (APV->Eta > 0) {
1299         if (APV->Thickness < 0.04) {
1300           (Charge_Vs_PathlengthTECP1[elepos])->Fill((*path)[i], Charge);
1301         }  // TEC plus, type 1
1302         else if (APV->Thickness > 0.04) {
1303           (Charge_Vs_PathlengthTECP2[elepos])->Fill((*path)[i], Charge);
1304         }  // TEC plus, type 2
1305       }
1306     }
1307 
1308   }  // END OF ON-CLUSTER LOOP
1309 }  //END OF processEvent()
1310 
1311 void SiStripGainFromCalibTree::algoAnalyzeTheTree() {
1312   for (unsigned int i = 0; i < VInputFiles.size(); i++) {
1313     printf("Openning file %3i/%3i --> %s\n", i + 1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str()));
1314     fflush(stdout);
1315     TFile* tfile = TFile::Open(VInputFiles[i].c_str());
1316     TString tree_path = TString::Format("gainCalibrationTree%s/tree", m_calibrationMode.c_str());
1317     TTree* tree = dynamic_cast<TTree*>(tfile->Get(tree_path.Data()));
1318 
1319     tree->SetBranchAddress((EventPrefix_ + "event" + EventSuffix_).c_str(), &eventnumber, nullptr);
1320     tree->SetBranchAddress((EventPrefix_ + "run" + EventSuffix_).c_str(), &runnumber, nullptr);
1321     tree->SetBranchAddress((EventPrefix_ + "TrigTech" + EventSuffix_).c_str(), &TrigTech, nullptr);
1322 
1323     tree->SetBranchAddress((TrackPrefix_ + "chi2ndof" + TrackSuffix_).c_str(), &trackchi2ndof, nullptr);
1324     tree->SetBranchAddress((TrackPrefix_ + "momentum" + TrackSuffix_).c_str(), &trackp, nullptr);
1325     tree->SetBranchAddress((TrackPrefix_ + "pt" + TrackSuffix_).c_str(), &trackpt, nullptr);
1326     tree->SetBranchAddress((TrackPrefix_ + "eta" + TrackSuffix_).c_str(), &tracketa, nullptr);
1327     tree->SetBranchAddress((TrackPrefix_ + "phi" + TrackSuffix_).c_str(), &trackphi, nullptr);
1328     tree->SetBranchAddress((TrackPrefix_ + "hitsvalid" + TrackSuffix_).c_str(), &trackhitsvalid, nullptr);
1329     tree->SetBranchAddress((TrackPrefix_ + "algo" + TrackSuffix_).c_str(), &trackalgo, nullptr);
1330 
1331     tree->SetBranchAddress((CalibPrefix_ + "trackindex" + CalibSuffix_).c_str(), &trackindex, nullptr);
1332     tree->SetBranchAddress((CalibPrefix_ + "rawid" + CalibSuffix_).c_str(), &rawid, nullptr);
1333     tree->SetBranchAddress((CalibPrefix_ + "localdirx" + CalibSuffix_).c_str(), &localdirx, nullptr);
1334     tree->SetBranchAddress((CalibPrefix_ + "localdiry" + CalibSuffix_).c_str(), &localdiry, nullptr);
1335     tree->SetBranchAddress((CalibPrefix_ + "localdirz" + CalibSuffix_).c_str(), &localdirz, nullptr);
1336     tree->SetBranchAddress((CalibPrefix_ + "firststrip" + CalibSuffix_).c_str(), &firststrip, nullptr);
1337     tree->SetBranchAddress((CalibPrefix_ + "nstrips" + CalibSuffix_).c_str(), &nstrips, nullptr);
1338     tree->SetBranchAddress((CalibPrefix_ + "saturation" + CalibSuffix_).c_str(), &saturation, nullptr);
1339     tree->SetBranchAddress((CalibPrefix_ + "overlapping" + CalibSuffix_).c_str(), &overlapping, nullptr);
1340     tree->SetBranchAddress((CalibPrefix_ + "farfromedge" + CalibSuffix_).c_str(), &farfromedge, nullptr);
1341     tree->SetBranchAddress((CalibPrefix_ + "charge" + CalibSuffix_).c_str(), &charge, nullptr);
1342     tree->SetBranchAddress((CalibPrefix_ + "path" + CalibSuffix_).c_str(), &path, nullptr);
1343     tree->SetBranchAddress((CalibPrefix_ + "chargeoverpath" + CalibSuffix_).c_str(), &chargeoverpath, nullptr);
1344     tree->SetBranchAddress((CalibPrefix_ + "amplitude" + CalibSuffix_).c_str(), &amplitude, nullptr);
1345     tree->SetBranchAddress((CalibPrefix_ + "gainused" + CalibSuffix_).c_str(), &gainused, nullptr);
1346     tree->SetBranchAddress((CalibPrefix_ + "gainusedTick" + CalibSuffix_).c_str(), &gainusedTick, nullptr);
1347 
1348     unsigned int nentries = tree->GetEntries();
1349     printf("Number of Events = %i + %i = %i\n", NEvent, nentries, (NEvent + nentries));
1350     printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
1351     printf("Looping on the Tree          :");
1352     int TreeStep = nentries / 50;
1353     if (TreeStep <= 1)
1354       TreeStep = 1;
1355     for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
1356       if (ientry % TreeStep == 0) {
1357         printf(".");
1358         fflush(stdout);
1359       }
1360       tree->GetEntry(ientry);
1361       processEvent();
1362     }
1363     printf("\n");  // END OF EVENT LOOP
1364   }
1365 }
1366 
1367 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
1368   unsigned int I = 0;
1369   TH1F* Proj = nullptr;
1370   double FitResults[6];
1371   double MPVmean = 300;
1372 
1373   int elepos = (AlgoMode == "PCL") ? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
1374 
1375   if (Charge_Vs_Index[elepos] == nullptr) {
1376     edm::LogError("SiStripGainFromCalibTree")
1377         << "Harvesting: could not execute algoComputeMPVandGain method because " << m_calibrationMode.c_str()
1378         << " statistics cannot be retrieved.\n"
1379         << "Please check if input contains " << m_calibrationMode.c_str() << " data." << std::endl;
1380     return;
1381   }
1382 
1383   TH2S* chvsidx = (Charge_Vs_Index[elepos])->getTH2S();
1384 
1385   printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
1386   printf("Fitting Charge Distribution  :");
1387   int TreeStep = APVsColl.size() / 50;
1388   for (auto it = APVsColl.begin(); it != APVsColl.end(); it++, I++) {
1389     if (I % TreeStep == 0) {
1390       printf(".");
1391       fflush(stdout);
1392     }
1393     stAPVGain* APV = it->second;
1394     if (APV->Bin < 0)
1395       APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
1396 
1397     if (APV->isMasked) {
1398       APV->Gain = APV->PreviousGain;
1399       MASKED++;
1400       continue;
1401     }
1402 
1403     Proj = (TH1F*)(chvsidx->ProjectionY(
1404         "", chvsidx->GetXaxis()->FindBin(APV->Index), chvsidx->GetXaxis()->FindBin(APV->Index), "e"));
1405     if (!Proj)
1406       continue;
1407 
1408     if (CalibrationLevel == 0) {
1409     } else if (CalibrationLevel == 1) {
1410       int SecondAPVId = APV->APVId;
1411       if (SecondAPVId % 2 == 0) {
1412         SecondAPVId = SecondAPVId + 1;
1413       } else {
1414         SecondAPVId = SecondAPVId - 1;
1415       }
1416       stAPVGain* APV2 = APVsColl[(APV->DetId << 4) | SecondAPVId];
1417       if (APV2->Bin < 0)
1418         APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1419       TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e"));
1420       if (Proj2) {
1421         Proj->Add(Proj2, 1);
1422         delete Proj2;
1423       }
1424     } else if (CalibrationLevel == 2) {
1425       for (unsigned int i = 0; i < 16; i++) {  //loop up to 6APV for Strip and up to 16 for Pixels
1426         auto tmpit = APVsColl.find((APV->DetId << 4) | i);
1427         if (tmpit == APVsColl.end())
1428           continue;
1429         stAPVGain* APV2 = tmpit->second;
1430         if (APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)
1431           continue;
1432         if (APV2->Bin < 0)
1433           APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1434         TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e"));
1435         if (Proj2) {
1436           Proj->Add(Proj2, 1);
1437           delete Proj2;
1438         }
1439       }
1440     } else {
1441       CalibrationLevel = 0;
1442       printf("Unknown Calibration Level, will assume %i\n", CalibrationLevel);
1443     }
1444 
1445     getPeakOfLandau(Proj, FitResults);
1446     APV->FitMPV = FitResults[0];
1447     APV->FitMPVErr = FitResults[1];
1448     APV->FitWidth = FitResults[2];
1449     APV->FitWidthErr = FitResults[3];
1450     APV->FitChi2 = FitResults[4];
1451     APV->FitNorm = FitResults[5];
1452     APV->NEntries = Proj->GetEntries();
1453 
1454     if (IsGoodLandauFit(FitResults)) {
1455       APV->Gain = APV->FitMPV / MPVmean;
1456       if (APV->SubDet > 2)
1457         GOOD++;
1458     } else {
1459       APV->Gain = APV->PreviousGain;
1460       if (APV->SubDet > 2)
1461         BAD++;
1462     }
1463     if (APV->Gain <= 0)
1464       APV->Gain = 1;
1465 
1466     //printf("%5i/%5i:  %6i - %1i  %5E Entries --> MPV = %f +- %f\n",I,APVsColl.size(),APV->DetId, APV->APVId, Proj->GetEntries(), FitResults[0], FitResults[1]);fflush(stdout);
1467     delete Proj;
1468   }
1469   printf("\n");
1470 }
1471 
1472 void SiStripGainFromCalibTree::qualityMonitor() {
1473   int elepos = (AlgoMode == "PCL") ? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
1474 
1475   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
1476     stAPVGain* APV = APVsCollOrdered[a];
1477     if (APV == nullptr)
1478       continue;
1479 
1480     unsigned int Index = APV->Index;
1481     unsigned int SubDet = APV->SubDet;
1482     unsigned int DetId = APV->DetId;
1483     float z = APV->z;
1484     float Eta = APV->Eta;
1485     float R = APV->R;
1486     float Phi = APV->Phi;
1487     float Thickness = APV->Thickness;
1488     double FitMPV = APV->FitMPV;
1489     double FitMPVErr = APV->FitMPVErr;
1490     double Gain = APV->Gain;
1491     double NEntries = APV->NEntries;
1492     double PreviousGain = APV->PreviousGain;
1493 
1494     if (SubDet < 3)
1495       continue;  // avoid to loop over Pixel det id
1496 
1497     if (Gain != 1.) {
1498       std::vector<MonitorElement*> charge_histos = APVGain::FetchMonitor(newCharge, DetId, tTopo_);
1499       TH2S* chvsidx = (Charge_Vs_Index[elepos])->getTH2S();
1500       int bin = chvsidx->GetXaxis()->FindBin(Index);
1501       TH1D* Proj = chvsidx->ProjectionY("proj", bin, bin);
1502       for (int binId = 0; binId < Proj->GetXaxis()->GetNbins(); binId++) {
1503         double new_charge = Proj->GetXaxis()->GetBinCenter(binId) / Gain;
1504         if (Proj->GetBinContent(binId) != 0.) {
1505           for (unsigned int h = 0; h < charge_histos.size(); h++) {
1506             TH1D* chisto = (charge_histos[h])->getTH1D();
1507             for (int e = 0; e < Proj->GetBinContent(binId); e++)
1508               chisto->Fill(new_charge);
1509           }
1510         }
1511       }
1512     }
1513 
1514     if (FitMPV <= 0.) {  // No fit of MPV
1515       if (APV->isMasked)
1516         NoMPVmasked->Fill(z, R);
1517       else
1518         NoMPVfit->Fill(z, R);
1519 
1520     } else {  // Fit of MPV
1521       if (FitMPV > 0.)
1522         Gains->Fill(Gain);
1523 
1524       MPVs->Fill(FitMPV);
1525       if (Thickness < 0.04)
1526         MPVs320->Fill(FitMPV);
1527       if (Thickness > 0.04)
1528         MPVs500->Fill(FitMPV);
1529 
1530       MPVError->Fill(FitMPVErr);
1531       MPVErrorVsMPV->Fill(FitMPV, FitMPVErr);
1532       MPVErrorVsEta->Fill(Eta, FitMPVErr);
1533       MPVErrorVsPhi->Fill(Phi, FitMPVErr);
1534       MPVErrorVsN->Fill(NEntries, FitMPVErr);
1535 
1536       if (SubDet == 3) {
1537         MPV_Vs_EtaTIB->Fill(Eta, FitMPV);
1538         MPV_Vs_PhiTIB->Fill(Phi, FitMPV);
1539         MPVsTIB->Fill(FitMPV);
1540 
1541       } else if (SubDet == 4) {
1542         MPV_Vs_EtaTID->Fill(Eta, FitMPV);
1543         MPV_Vs_PhiTID->Fill(Phi, FitMPV);
1544         MPVsTID->Fill(FitMPV);
1545         if (Eta < 0.)
1546           MPVsTIDM->Fill(FitMPV);
1547         if (Eta > 0.)
1548           MPVsTIDP->Fill(FitMPV);
1549 
1550       } else if (SubDet == 5) {
1551         MPV_Vs_EtaTOB->Fill(Eta, FitMPV);
1552         MPV_Vs_PhiTOB->Fill(Phi, FitMPV);
1553         MPVsTOB->Fill(FitMPV);
1554 
1555       } else if (SubDet == 6) {
1556         MPV_Vs_EtaTEC->Fill(Eta, FitMPV);
1557         MPV_Vs_PhiTEC->Fill(Phi, FitMPV);
1558         MPVsTEC->Fill(FitMPV);
1559         if (Eta < 0.)
1560           MPVsTECM->Fill(FitMPV);
1561         if (Eta > 0.)
1562           MPVsTECP->Fill(FitMPV);
1563         if (Thickness < 0.04) {
1564           MPV_Vs_EtaTECthin->Fill(Eta, FitMPV);
1565           MPV_Vs_PhiTECthin->Fill(Phi, FitMPV);
1566           MPVsTECthin->Fill(FitMPV);
1567           if (Eta > 0.)
1568             MPVsTECP1->Fill(FitMPV);
1569           if (Eta < 0.)
1570             MPVsTECM1->Fill(FitMPV);
1571         }
1572         if (Thickness > 0.04) {
1573           MPV_Vs_EtaTECthick->Fill(Eta, FitMPV);
1574           MPV_Vs_PhiTECthick->Fill(Phi, FitMPV);
1575           MPVsTECthick->Fill(FitMPV);
1576           if (Eta > 0.)
1577             MPVsTECP2->Fill(FitMPV);
1578           if (Eta < 0.)
1579             MPVsTECM2->Fill(FitMPV);
1580         }
1581       }
1582     }
1583 
1584     if (SubDet == 3 && PreviousGain != 0.)
1585       DiffWRTPrevGainTIB->Fill(Gain / PreviousGain);
1586     else if (SubDet == 4 && PreviousGain != 0.)
1587       DiffWRTPrevGainTID->Fill(Gain / PreviousGain);
1588     else if (SubDet == 5 && PreviousGain != 0.)
1589       DiffWRTPrevGainTOB->Fill(Gain / PreviousGain);
1590     else if (SubDet == 6 && PreviousGain != 0.)
1591       DiffWRTPrevGainTEC->Fill(Gain / PreviousGain);
1592 
1593     if (SubDet == 3)
1594       GainVsPrevGainTIB->Fill(PreviousGain, Gain);
1595     else if (SubDet == 4)
1596       GainVsPrevGainTID->Fill(PreviousGain, Gain);
1597     else if (SubDet == 5)
1598       GainVsPrevGainTOB->Fill(PreviousGain, Gain);
1599     else if (SubDet == 6)
1600       GainVsPrevGainTEC->Fill(PreviousGain, Gain);
1601   }
1602 }
1603 
1604 void SiStripGainFromCalibTree::storeOnTree(TFileService* tfs) {
1605   unsigned int tree_Index;
1606   unsigned int tree_Bin;
1607   unsigned int tree_DetId;
1608   unsigned char tree_APVId;
1609   unsigned char tree_SubDet;
1610   float tree_x;
1611   float tree_y;
1612   float tree_z;
1613   float tree_Eta;
1614   float tree_R;
1615   float tree_Phi;
1616   float tree_Thickness;
1617   float tree_FitMPV;
1618   float tree_FitMPVErr;
1619   float tree_FitWidth;
1620   float tree_FitWidthErr;
1621   float tree_FitChi2NDF;
1622   float tree_FitNorm;
1623   double tree_Gain;
1624   double tree_PrevGain;
1625   double tree_PrevGainTick;
1626   double tree_NEntries;
1627   bool tree_isMasked;
1628 
1629   TTree* MyTree;
1630   MyTree = tfs->make<TTree>("APVGain", "APVGain");
1631   MyTree->Branch("Index", &tree_Index, "Index/i");
1632   MyTree->Branch("Bin", &tree_Bin, "Bin/i");
1633   MyTree->Branch("DetId", &tree_DetId, "DetId/i");
1634   MyTree->Branch("APVId", &tree_APVId, "APVId/b");
1635   MyTree->Branch("SubDet", &tree_SubDet, "SubDet/b");
1636   MyTree->Branch("x", &tree_x, "x/F");
1637   MyTree->Branch("y", &tree_y, "y/F");
1638   MyTree->Branch("z", &tree_z, "z/F");
1639   MyTree->Branch("Eta", &tree_Eta, "Eta/F");
1640   MyTree->Branch("R", &tree_R, "R/F");
1641   MyTree->Branch("Phi", &tree_Phi, "Phi/F");
1642   MyTree->Branch("Thickness", &tree_Thickness, "Thickness/F");
1643   MyTree->Branch("FitMPV", &tree_FitMPV, "FitMPV/F");
1644   MyTree->Branch("FitMPVErr", &tree_FitMPVErr, "FitMPVErr/F");
1645   MyTree->Branch("FitWidth", &tree_FitWidth, "FitWidth/F");
1646   MyTree->Branch("FitWidthErr", &tree_FitWidthErr, "FitWidthErr/F");
1647   MyTree->Branch("FitChi2NDF", &tree_FitChi2NDF, "FitChi2NDF/F");
1648   MyTree->Branch("FitNorm", &tree_FitNorm, "FitNorm/F");
1649   MyTree->Branch("Gain", &tree_Gain, "Gain/D");
1650   MyTree->Branch("PrevGain", &tree_PrevGain, "PrevGain/D");
1651   MyTree->Branch("PrevGainTick", &tree_PrevGainTick, "PrevGainTick/D");
1652   MyTree->Branch("NEntries", &tree_NEntries, "NEntries/D");
1653   MyTree->Branch("isMasked", &tree_isMasked, "isMasked/O");
1654 
1655   FILE* Gains = stdout;
1656   fprintf(Gains, "NEvents   = %i\n", NEvent);
1657   fprintf(Gains, "NTracks   = %i\n", NTrack);
1658   //fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1659   fprintf(Gains, "NClustersStrip = %i\n", NClusterStrip);
1660   //fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1661   fprintf(Gains, "Number of Strip APVs = %lu\n", static_cast<unsigned long>(NStripAPVs));
1662   fprintf(Gains,
1663           "GoodFits = %i BadFits = %i ratio = %f%%   (MASKED=%i)\n",
1664           GOOD,
1665           BAD,
1666           (100.0 * GOOD) / (GOOD + BAD),
1667           MASKED);
1668 
1669   Gains = fopen(OutputGains.c_str(), "w");
1670   fprintf(Gains, "NEvents   = %i\n", NEvent);
1671   fprintf(Gains, "NTracks   = %i\n", NTrack);
1672   //fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1673   fprintf(Gains, "NClustersStrip = %i\n", NClusterStrip);
1674   fprintf(Gains, "Number of Strip APVs = %lu\n", static_cast<unsigned long>(NStripAPVs));
1675   //fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1676   fprintf(Gains,
1677           "GoodFits = %i BadFits = %i ratio = %f%%   (MASKED=%i)\n",
1678           GOOD,
1679           BAD,
1680           (100.0 * GOOD) / (GOOD + BAD),
1681           MASKED);
1682 
1683   int elepos = statCollectionFromMode(m_calibrationMode.c_str());
1684 
1685   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
1686     stAPVGain* APV = APVsCollOrdered[a];
1687     if (APV == nullptr)
1688       continue;
1689     //     printf(      "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
1690     fprintf(Gains,
1691             "%i | %i | PreviousGain = %7.5f(tick) x %7.5f(particle) NewGain (particle) = %7.5f (#clusters=%8.0f)\n",
1692             APV->DetId,
1693             APV->APVId,
1694             APV->PreviousGainTick,
1695             APV->PreviousGain,
1696             APV->Gain,
1697             APV->NEntries);
1698 
1699     tree_Index = APV->Index;
1700     tree_Bin = (Charge_Vs_Index[elepos])->getTH2S()->GetXaxis()->FindBin(APV->Index);
1701     tree_DetId = APV->DetId;
1702     tree_APVId = APV->APVId;
1703     tree_SubDet = APV->SubDet;
1704     tree_x = APV->x;
1705     tree_y = APV->y;
1706     tree_z = APV->z;
1707     tree_Eta = APV->Eta;
1708     tree_R = APV->R;
1709     tree_Phi = APV->Phi;
1710     tree_Thickness = APV->Thickness;
1711     tree_FitMPV = APV->FitMPV;
1712     tree_FitMPVErr = APV->FitMPVErr;
1713     tree_FitWidth = APV->FitWidth;
1714     tree_FitWidthErr = APV->FitWidthErr;
1715     tree_FitChi2NDF = APV->FitChi2;
1716     tree_FitNorm = APV->FitNorm;
1717     tree_Gain = APV->Gain;
1718     tree_PrevGain = APV->PreviousGain;
1719     tree_PrevGainTick = APV->PreviousGainTick;
1720     tree_NEntries = APV->NEntries;
1721     tree_isMasked = APV->isMasked;
1722 
1723     if (tree_DetId == 402673324) {
1724       printf("%i | %i : %f --> %f  (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
1725     }
1726 
1727     MyTree->Fill();
1728   }
1729   if (Gains)
1730     fclose(Gains);
1731 }
1732 
1733 bool SiStripGainFromCalibTree::produceTagFilter() {
1734   // The goal of this function is to check wether or not there is enough statistics
1735   // to produce a meaningful tag for the DB
1736   int elepos = (AlgoMode == "PCL") ? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
1737   if (Charge_Vs_Index[elepos] == nullptr) {
1738     edm::LogError("SiStripGainFromCalibTree")
1739         << "produceTagFilter -> Return false: could not retrieve the " << m_calibrationMode.c_str() << " statistics.\n"
1740         << "Please check if input contains " << m_calibrationMode.c_str() << " data." << std::endl;
1741     return false;
1742   }
1743 
1744   float integral = (Charge_Vs_Index[elepos])->getTH2S()->Integral();
1745   if ((Charge_Vs_Index[elepos])->getTH2S()->Integral(0, NStripAPVs + 1, 0, 99999) < tagCondition_NClusters) {
1746     edm::LogWarning("SiStripGainFromCalibTree")
1747         << "calibrationMode  -> " << m_calibrationMode << "\n"
1748         << "produceTagFilter -> Return false: Statistics is too low : " << integral << endl;
1749     return false;
1750   }
1751   if ((1.0 * GOOD) / (GOOD + BAD) < tagCondition_GoodFrac) {
1752     edm::LogWarning("SiStripGainFromCalibTree")
1753         << "calibrationMode  -> " << m_calibrationMode << "\n"
1754         << "produceTagFilter ->  Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD + BAD) << endl;
1755     return false;
1756   }
1757   return true;
1758 }
1759 
1760 std::unique_ptr<SiStripApvGain> SiStripGainFromCalibTree::getNewObject() {
1761   auto obj = std::make_unique<SiStripApvGain>();
1762   if (!m_harvestingMode)
1763     return obj;
1764 
1765   if (!produceTagFilter()) {
1766     edm::LogWarning("SiStripGainFromCalibTree")
1767         << "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
1768     setDoStore(false);
1769     return obj;
1770   }
1771 
1772   std::vector<float>* theSiStripVector = nullptr;
1773   unsigned int PreviousDetId = 0;
1774   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
1775     stAPVGain* APV = APVsCollOrdered[a];
1776     if (APV == nullptr) {
1777       printf("Bug\n");
1778       continue;
1779     }
1780     if (APV->SubDet <= 2)
1781       continue;
1782     if (APV->DetId != PreviousDetId) {
1783       if (theSiStripVector != nullptr) {
1784         SiStripApvGain::Range range(theSiStripVector->begin(), theSiStripVector->end());
1785         if (!obj->put(PreviousDetId, range))
1786           printf("Bug to put detId = %i\n", PreviousDetId);
1787       }
1788       theSiStripVector = new std::vector<float>;
1789       PreviousDetId = APV->DetId;
1790     }
1791     if (theSiStripVector != nullptr) {
1792       theSiStripVector->push_back(APV->Gain);
1793     }
1794   }
1795   if (theSiStripVector != nullptr) {
1796     SiStripApvGain::Range range(theSiStripVector->begin(), theSiStripVector->end());
1797     if (!obj->put(PreviousDetId, range))
1798       printf("Bug to put detId = %i\n", PreviousDetId);
1799   }
1800 
1801   if (theSiStripVector != nullptr)
1802     delete theSiStripVector;
1803 
1804   return obj;
1805 }
1806 
1807 SiStripGainFromCalibTree::~SiStripGainFromCalibTree() {
1808   APVsColl.clear();
1809   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
1810     stAPVGain* APV = APVsCollOrdered[a];
1811     if (APV != nullptr)
1812       delete APV;
1813   }
1814   APVsCollOrdered.clear();
1815 }
1816 
1817 void SiStripGainFromCalibTree::MakeCalibrationMap() {
1818   if (!useCalibration)
1819     return;
1820 
1821   TChain* t1 = new TChain("SiStripCalib/APVGain");
1822   t1->Add(m_calibrationPath.c_str());
1823 
1824   unsigned int tree_DetId;
1825   unsigned char tree_APVId;
1826   double tree_Gain;
1827 
1828   t1->SetBranchAddress("DetId", &tree_DetId);
1829   t1->SetBranchAddress("APVId", &tree_APVId);
1830   t1->SetBranchAddress("Gain", &tree_Gain);
1831 
1832   for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
1833     t1->GetEntry(ientry);
1834     stAPVGain* APV = APVsColl[(tree_DetId << 4) | (unsigned int)tree_APVId];
1835     APV->CalibGain = tree_Gain;
1836   }
1837 
1838   delete t1;
1839 }
1840 
1841 void SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1842   // in AlCaHarvesting mode we just need to run the logic in the endJob step
1843   if (m_harvestingMode)
1844     return;
1845 
1846   if (AlgoMode == "CalibTree")
1847     return;
1848 
1849   eventnumber = iEvent.id().event();
1850   runnumber = iEvent.id().run();
1851   auto handle01 = connect(TrigTech, TrigTech_token_, iEvent);
1852   auto handle02 = connect(trackchi2ndof, trackchi2ndof_token_, iEvent);
1853   auto handle03 = connect(trackp, trackp_token_, iEvent);
1854   auto handle04 = connect(trackpt, trackpt_token_, iEvent);
1855   auto handle05 = connect(tracketa, tracketa_token_, iEvent);
1856   auto handle06 = connect(trackphi, trackphi_token_, iEvent);
1857   auto handle07 = connect(trackhitsvalid, trackhitsvalid_token_, iEvent);
1858   auto handle08 = connect(trackindex, trackindex_token_, iEvent);
1859   auto handle09 = connect(rawid, rawid_token_, iEvent);
1860   auto handle11 = connect(localdirx, localdirx_token_, iEvent);
1861   auto handle12 = connect(localdiry, localdiry_token_, iEvent);
1862   auto handle13 = connect(localdirz, localdirz_token_, iEvent);
1863   auto handle14 = connect(firststrip, firststrip_token_, iEvent);
1864   auto handle15 = connect(nstrips, nstrips_token_, iEvent);
1865   auto handle16 = connect(saturation, saturation_token_, iEvent);
1866   auto handle17 = connect(overlapping, overlapping_token_, iEvent);
1867   auto handle18 = connect(farfromedge, farfromedge_token_, iEvent);
1868   auto handle19 = connect(charge, charge_token_, iEvent);
1869   auto handle21 = connect(path, path_token_, iEvent);
1870   auto handle22 = connect(chargeoverpath, chargeoverpath_token_, iEvent);
1871   auto handle23 = connect(amplitude, amplitude_token_, iEvent);
1872   auto handle24 = connect(gainused, gainused_token_, iEvent);
1873   auto handle25 = connect(gainusedTick, gainusedTick_token_, iEvent);
1874   auto handle26 = connect(trackalgo, trackalgo_token_, iEvent);
1875 
1876   processEvent();
1877 }
1878 
1879 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);