Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-25 05:34:36

0001 // system include files
0002 #include <memory>
0003 
0004 // CMSSW includes
0005 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
0006 #include "DataFormats/DetId/interface/DetId.h"
0007 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
0008 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0012 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0015 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0016 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0017 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0018 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0019 
0020 // user include files
0021 #include "CalibTracker/SiStripChannelGain/interface/SiStripGainsPCLHarvester.h"
0022 #include "CalibTracker/SiStripChannelGain/interface/APVGainHelpers.h"
0023 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0024 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include <iostream>
0027 #include <sstream>
0028 
0029 // ROOT includes
0030 #include <TROOT.h>
0031 #include <TTree.h>
0032 
0033 //********************************************************************************//
0034 SiStripGainsPCLHarvester::SiStripGainsPCLHarvester(const edm::ParameterSet& ps)
0035     : doStoreOnDB(false), GOOD(0), BAD(0), MASKED(0), NStripAPVs(0), NPixelDets(0) {
0036   m_Record = ps.getUntrackedParameter<std::string>("Record", "SiStripApvGainRcd");
0037   CalibrationLevel = ps.getUntrackedParameter<int>("CalibrationLevel", 0);
0038   MinNrEntries = ps.getUntrackedParameter<double>("minNrEntries", 20);
0039   m_DQMdir = ps.getUntrackedParameter<std::string>("DQMdir", "AlCaReco/SiStripGains");
0040   m_calibrationMode = ps.getUntrackedParameter<std::string>("calibrationMode", "StdBunch");
0041   tagCondition_NClusters = ps.getUntrackedParameter<double>("NClustersForTagProd", 2E8);
0042   tagCondition_GoodFrac = ps.getUntrackedParameter<double>("GoodFracForTagProd", 0.95);
0043   doChargeMonitorPerPlane = ps.getUntrackedParameter<bool>("doChargeMonitorPerPlane", false);
0044   VChargeHisto = ps.getUntrackedParameter<std::vector<std::string>>("ChargeHisto");
0045   fit_gaussianConvolution_ = ps.getUntrackedParameter<bool>("FitGaussianConvolution", false);
0046   fit_gaussianConvolutionTOBL56_ = ps.getUntrackedParameter<bool>("FitGaussianConvolutionTOBL5L6", false);
0047   fit_dataDrivenRange_ = ps.getUntrackedParameter<bool>("FitDataDrivenRange", false);
0048   storeGainsTree_ = ps.getUntrackedParameter<bool>("StoreGainsTree", false);
0049 
0050   //Set the monitoring element tag and store
0051   dqm_tag_.reserve(7);
0052   dqm_tag_.clear();
0053   dqm_tag_.push_back("StdBunch");    // statistic collection from Standard Collision Bunch @ 3.8 T
0054   dqm_tag_.push_back("StdBunch0T");  // statistic collection from Standard Collision Bunch @ 0 T
0055   dqm_tag_.push_back("AagBunch");    // statistic collection from First Collision After Abort Gap @ 3.8 T
0056   dqm_tag_.push_back("AagBunch0T");  // statistic collection from First Collision After Abort Gap @ 0 T
0057   dqm_tag_.push_back("IsoMuon");     // statistic collection from Isolated Muon @ 3.8 T
0058   dqm_tag_.push_back("IsoMuon0T");   // statistic collection from Isolated Muon @ 0 T
0059   dqm_tag_.push_back("Harvest");     // statistic collection: Harvest
0060 
0061   tTopoToken_ = esConsumes<edm::Transition::EndRun>();
0062   tkGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0063   gainToken_ = esConsumes<edm::Transition::BeginRun>();
0064   qualityToken_ = esConsumes<edm::Transition::BeginRun>();
0065 }
0066 
0067 //********************************************************************************//
0068 // ------------ method called for each event  ------------
0069 void SiStripGainsPCLHarvester::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0070   using namespace edm;
0071   static constexpr float defaultGainTick = 690. / 640.;
0072 
0073   this->checkBookAPVColls(iSetup);  // check whether APV colls are booked and do so if not yet done
0074 
0075   const auto gainHandle = iSetup.getHandle(gainToken_);
0076   if (!gainHandle.isValid()) {
0077     edm::LogError("SiStripGainPCLHarvester") << "gainHandle is not valid\n";
0078     exit(0);
0079   }
0080 
0081   const auto& stripQuality = iSetup.getData(qualityToken_);
0082 
0083   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0084     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0085 
0086     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0087       continue;
0088 
0089     APV->isMasked = stripQuality.IsApvBad(APV->DetId, APV->APVId);
0090 
0091     if (gainHandle->getNumberOfTags() != 2) {
0092       edm::LogError("SiStripGainPCLHarvester") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
0093       fflush(stdout);
0094       exit(0);
0095     };
0096     float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
0097     if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
0098       edm::LogWarning("SiStripGainPCLHarvester") << "WARNING: ParticleGain in the global tag changed\n";
0099     APV->PreviousGain = newPreviousGain;
0100 
0101     float newPreviousGainTick =
0102         APV->isMasked ? defaultGainTick : gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
0103     if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
0104       edm::LogWarning("SiStripGainPCLHarvester")
0105           << "WARNING: TickMarkGain in the global tag changed\n"
0106           << std::endl
0107           << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
0108           << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
0109           << std::endl;
0110     }
0111     APV->PreviousGainTick = newPreviousGainTick;
0112   }
0113 }
0114 
0115 //********************************************************************************//
0116 void SiStripGainsPCLHarvester::dqmEndJob(DQMStore::IBooker& ibooker_, DQMStore::IGetter& igetter_) {
0117   edm::LogInfo("SiStripGainsPCLHarvester") << "Starting harvesting statistics" << std::endl;
0118 
0119   std::string DQM_dir = m_DQMdir;
0120 
0121   std::string stag = *(std::find(dqm_tag_.begin(), dqm_tag_.end(), m_calibrationMode));
0122   if (!stag.empty() && stag[0] != '_')
0123     stag.insert(0, 1, '_');
0124 
0125   std::string cvi = DQM_dir + std::string("/Charge_Vs_Index") + stag;
0126 
0127   MonitorElement* Charge_Vs_Index = igetter_.get(cvi);
0128 
0129   if (Charge_Vs_Index == nullptr) {
0130     edm::LogError("SiStripGainsPCLHarvester")
0131         << "Harvesting: could not retrieve " << cvi.c_str() << ", statistics will not be summed!" << std::endl;
0132   } else {
0133     edm::LogInfo("SiStripGainsPCLHarvester")
0134         << "Harvesting " << (Charge_Vs_Index)->getTH2S()->GetEntries() << " more clusters" << std::endl;
0135   }
0136 
0137   algoComputeMPVandGain(Charge_Vs_Index);
0138   std::unique_ptr<SiStripApvGain> theAPVGains = this->getNewObject(Charge_Vs_Index);
0139 
0140   // write out the APVGains record
0141   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0142 
0143   if (doStoreOnDB) {
0144     if (poolDbService.isAvailable())
0145       poolDbService->writeOneIOV(*theAPVGains, poolDbService->currentTime(), m_Record);
0146     else
0147       throw std::runtime_error("PoolDBService required.");
0148   } else {
0149     edm::LogInfo("SiStripGainsPCLHarvester") << "Will not produce payload!" << std::endl;
0150   }
0151 
0152   //Collect the statistics for monitoring and validation
0153   gainQualityMonitor(ibooker_, Charge_Vs_Index);
0154 
0155   if (storeGainsTree_) {
0156     if (Charge_Vs_Index != nullptr) {
0157       storeGainsTree(Charge_Vs_Index->getTH2S()->GetXaxis());
0158     } else {
0159       edm::LogError("SiStripGainsPCLHarvester")
0160           << "Harvesting: could not retrieve " << cvi.c_str() << "Tree won't be stored" << std::endl;
0161     }
0162   }
0163 }
0164 
0165 //********************************************************************************//
0166 void SiStripGainsPCLHarvester::gainQualityMonitor(DQMStore::IBooker& ibooker_,
0167                                                   const MonitorElement* Charge_Vs_Index) const {
0168   ibooker_.setCurrentFolder("AlCaReco/SiStripGainsHarvesting/");
0169 
0170   std::vector<APVGain::APVmon> new_charge_histos;
0171   std::vector<std::pair<std::string, std::string>> cnames =
0172       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "newG2");
0173   for (unsigned int i = 0; i < cnames.size(); i++) {
0174     MonitorElement* monitor = ibooker_.book1DD((cnames[i]).first, (cnames[i]).second.c_str(), 100, 0., 1000.);
0175     int thick = APVGain::thickness((cnames[i]).first);
0176     int id = APVGain::subdetectorId((cnames[i]).first);
0177     int side = APVGain::subdetectorSide((cnames[i]).first);
0178     int plane = APVGain::subdetectorPlane((cnames[i]).first);
0179     new_charge_histos.push_back(APVGain::APVmon(thick, id, side, plane, monitor));
0180   }
0181 
0182   int MPVbin = 300;
0183   float MPVmin = 0.;
0184   float MPVmax = 600.;
0185 
0186   MonitorElement* MPV_Vs_EtaTIB =
0187       ibooker_.book2DD("MPVvsEtaTIB", "MPV vs Eta TIB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0188   MonitorElement* MPV_Vs_EtaTID =
0189       ibooker_.book2DD("MPVvsEtaTID", "MPV vs Eta TID", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0190   MonitorElement* MPV_Vs_EtaTOB =
0191       ibooker_.book2DD("MPVvsEtaTOB", "MPV vs Eta TOB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0192   MonitorElement* MPV_Vs_EtaTEC =
0193       ibooker_.book2DD("MPVvsEtaTEC", "MPV vs Eta TEC", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0194   MonitorElement* MPV_Vs_EtaTECthin =
0195       ibooker_.book2DD("MPVvsEtaTEC1", "MPV vs Eta TEC-thin", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0196   MonitorElement* MPV_Vs_EtaTECthick =
0197       ibooker_.book2DD("MPVvsEtaTEC2", "MPV vs Eta TEC-thick", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
0198 
0199   MonitorElement* MPV_Vs_PhiTIB =
0200       ibooker_.book2DD("MPVvsPhiTIB", "MPV vs Phi TIB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0201   MonitorElement* MPV_Vs_PhiTID =
0202       ibooker_.book2DD("MPVvsPhiTID", "MPV vs Phi TID", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0203   MonitorElement* MPV_Vs_PhiTOB =
0204       ibooker_.book2DD("MPVvsPhiTOB", "MPV vs Phi TOB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0205   MonitorElement* MPV_Vs_PhiTEC =
0206       ibooker_.book2DD("MPVvsPhiTEC", "MPV vs Phi TEC", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0207   MonitorElement* MPV_Vs_PhiTECthin =
0208       ibooker_.book2DD("MPVvsPhiTEC1", "MPV vs Phi TEC-thin ", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0209   MonitorElement* MPV_Vs_PhiTECthick =
0210       ibooker_.book2DD("MPVvsPhiTEC2", "MPV vs Phi TEC-thick", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
0211 
0212   MonitorElement* NoMPVfit = ibooker_.book2DD("NoMPVfit", "Modules with bad Landau Fit", 350, -350, 350, 240, 0, 120);
0213   MonitorElement* NoMPVmasked = ibooker_.book2DD("NoMPVmasked", "Masked Modules", 350, -350, 350, 240, 0, 120);
0214 
0215   MonitorElement* Gains = ibooker_.book1DD("Gains", "Gains", 300, 0, 2);
0216   MonitorElement* MPVs = ibooker_.book1DD("MPVs", "MPVs", MPVbin, MPVmin, MPVmax);
0217   MonitorElement* MPVs320 = ibooker_.book1DD("MPV_320", "MPV 320 thickness", MPVbin, MPVmin, MPVmax);
0218   MonitorElement* MPVs500 = ibooker_.book1DD("MPV_500", "MPV 500 thickness", MPVbin, MPVmin, MPVmax);
0219   MonitorElement* MPVsTIB = ibooker_.book1DD("MPV_TIB", "MPV TIB", MPVbin, MPVmin, MPVmax);
0220   MonitorElement* MPVsTID = ibooker_.book1DD("MPV_TID", "MPV TID", MPVbin, MPVmin, MPVmax);
0221   MonitorElement* MPVsTIDP = ibooker_.book1DD("MPV_TIDP", "MPV TIDP", MPVbin, MPVmin, MPVmax);
0222   MonitorElement* MPVsTIDM = ibooker_.book1DD("MPV_TIDM", "MPV TIDM", MPVbin, MPVmin, MPVmax);
0223   MonitorElement* MPVsTOB = ibooker_.book1DD("MPV_TOB", "MPV TOB", MPVbin, MPVmin, MPVmax);
0224   MonitorElement* MPVsTEC = ibooker_.book1DD("MPV_TEC", "MPV TEC", MPVbin, MPVmin, MPVmax);
0225   MonitorElement* MPVsTECP = ibooker_.book1DD("MPV_TECP", "MPV TECP", MPVbin, MPVmin, MPVmax);
0226   MonitorElement* MPVsTECM = ibooker_.book1DD("MPV_TECM", "MPV TECM", MPVbin, MPVmin, MPVmax);
0227   MonitorElement* MPVsTECthin = ibooker_.book1DD("MPV_TEC1", "MPV TEC thin", MPVbin, MPVmin, MPVmax);
0228   MonitorElement* MPVsTECthick = ibooker_.book1DD("MPV_TEC2", "MPV TEC thick", MPVbin, MPVmin, MPVmax);
0229   MonitorElement* MPVsTECP1 = ibooker_.book1DD("MPV_TECP1", "MPV TECP thin ", MPVbin, MPVmin, MPVmax);
0230   MonitorElement* MPVsTECP2 = ibooker_.book1DD("MPV_TECP2", "MPV TECP thick", MPVbin, MPVmin, MPVmax);
0231   MonitorElement* MPVsTECM1 = ibooker_.book1DD("MPV_TECM1", "MPV TECM thin", MPVbin, MPVmin, MPVmax);
0232   MonitorElement* MPVsTECM2 = ibooker_.book1DD("MPV_TECM2", "MPV TECM thick", MPVbin, MPVmin, MPVmax);
0233 
0234   MonitorElement* MPVError = ibooker_.book1DD("MPVError", "MPV Error", 150, 0, 150);
0235   MonitorElement* MPVErrorVsMPV = ibooker_.book2DD("MPVErrorVsMPV", "MPV Error vs MPV", 300, 0, 600, 150, 0, 150);
0236   MonitorElement* MPVErrorVsEta = ibooker_.book2DD("MPVErrorVsEta", "MPV Error vs Eta", 50, -3.0, 3.0, 150, 0, 150);
0237   MonitorElement* MPVErrorVsPhi = ibooker_.book2DD("MPVErrorVsPhi", "MPV Error vs Phi", 50, -3.4, 3.4, 150, 0, 150);
0238   MonitorElement* MPVErrorVsN = ibooker_.book2DD("MPVErrorVsN", "MPV Error vs N", 500, 0, 1000, 150, 0, 150);
0239 
0240   MonitorElement* DiffWRTPrevGainTIB = ibooker_.book1DD("DiffWRTPrevGainTIB", "Diff w.r.t. PrevGain TIB", 250, 0, 2);
0241   MonitorElement* DiffWRTPrevGainTID = ibooker_.book1DD("DiffWRTPrevGainTID", "Diff w.r.t. PrevGain TID", 250, 0, 2);
0242   MonitorElement* DiffWRTPrevGainTOB = ibooker_.book1DD("DiffWRTPrevGainTOB", "Diff w.r.t. PrevGain TOB", 250, 0, 2);
0243   MonitorElement* DiffWRTPrevGainTEC = ibooker_.book1DD("DiffWRTPrevGainTEC", "Diff w.r.t. PrevGain TEC", 250, 0, 2);
0244 
0245   MonitorElement* GainVsPrevGainTIB =
0246       ibooker_.book2DD("GainVsPrevGainTIB", "Gain vs PrevGain TIB", 100, 0, 2, 100, 0, 2);
0247   MonitorElement* GainVsPrevGainTID =
0248       ibooker_.book2DD("GainVsPrevGainTID", "Gain vs PrevGain TID", 100, 0, 2, 100, 0, 2);
0249   MonitorElement* GainVsPrevGainTOB =
0250       ibooker_.book2DD("GainVsPrevGainTOB", "Gain vs PrevGain TOB", 100, 0, 2, 100, 0, 2);
0251   MonitorElement* GainVsPrevGainTEC =
0252       ibooker_.book2DD("GainVsPrevGainTEC", "Gain vs PrevGain TEC", 100, 0, 2, 100, 0, 2);
0253 
0254   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0255     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0256     if (APV == nullptr)
0257       continue;
0258 
0259     unsigned int Index = APV->Index;
0260     unsigned int SubDet = APV->SubDet;
0261     unsigned int DetId = APV->DetId;
0262     float z = APV->z;
0263     float Eta = APV->Eta;
0264     float R = APV->R;
0265     float Phi = APV->Phi;
0266     float Thickness = APV->Thickness;
0267     double FitMPV = APV->FitMPV;
0268     double FitMPVErr = APV->FitMPVErr;
0269     double Gain = APV->Gain;
0270     double NEntries = APV->NEntries;
0271     double PreviousGain = APV->PreviousGain;
0272 
0273     if (SubDet < 3)
0274       continue;  // avoid to loop over Pixel det id
0275 
0276     if (Gain != 1.) {
0277       std::vector<MonitorElement*> charge_histos = APVGain::FetchMonitor(new_charge_histos, DetId, tTopo_.get());
0278 
0279       if (!Charge_Vs_Index)
0280         continue;
0281       TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
0282       int bin = chvsidx->GetXaxis()->FindBin(Index);
0283       TH1D* Proj = chvsidx->ProjectionY("proj", bin, bin);
0284       for (int binId = 0; binId < Proj->GetXaxis()->GetNbins(); binId++) {
0285         double new_charge = Proj->GetXaxis()->GetBinCenter(binId) / Gain;
0286         if (Proj->GetBinContent(binId) != 0.) {
0287           for (unsigned int h = 0; h < charge_histos.size(); h++) {
0288             TH1D* chisto = (charge_histos[h])->getTH1D();
0289             for (int e = 0; e < Proj->GetBinContent(binId); e++)
0290               chisto->Fill(new_charge);
0291           }
0292         }
0293       }
0294     }
0295 
0296     if (FitMPV <= 0.) {  // No fit of MPV
0297       if (APV->isMasked)
0298         NoMPVmasked->Fill(z, R);
0299       else
0300         NoMPVfit->Fill(z, R);
0301 
0302     } else {  // Fit of MPV
0303       if (FitMPV > 0.)
0304         Gains->Fill(Gain);
0305 
0306       MPVs->Fill(FitMPV);
0307       if (Thickness < 0.04)
0308         MPVs320->Fill(FitMPV);
0309       if (Thickness > 0.04)
0310         MPVs500->Fill(FitMPV);
0311 
0312       MPVError->Fill(FitMPVErr);
0313       MPVErrorVsMPV->Fill(FitMPV, FitMPVErr);
0314       MPVErrorVsEta->Fill(Eta, FitMPVErr);
0315       MPVErrorVsPhi->Fill(Phi, FitMPVErr);
0316       MPVErrorVsN->Fill(NEntries, FitMPVErr);
0317 
0318       if (SubDet == 3) {
0319         MPV_Vs_EtaTIB->Fill(Eta, FitMPV);
0320         MPV_Vs_PhiTIB->Fill(Phi, FitMPV);
0321         MPVsTIB->Fill(FitMPV);
0322 
0323       } else if (SubDet == 4) {
0324         MPV_Vs_EtaTID->Fill(Eta, FitMPV);
0325         MPV_Vs_PhiTID->Fill(Phi, FitMPV);
0326         MPVsTID->Fill(FitMPV);
0327         if (Eta < 0.)
0328           MPVsTIDM->Fill(FitMPV);
0329         if (Eta > 0.)
0330           MPVsTIDP->Fill(FitMPV);
0331 
0332       } else if (SubDet == 5) {
0333         MPV_Vs_EtaTOB->Fill(Eta, FitMPV);
0334         MPV_Vs_PhiTOB->Fill(Phi, FitMPV);
0335         MPVsTOB->Fill(FitMPV);
0336 
0337       } else if (SubDet == 6) {
0338         MPV_Vs_EtaTEC->Fill(Eta, FitMPV);
0339         MPV_Vs_PhiTEC->Fill(Phi, FitMPV);
0340         MPVsTEC->Fill(FitMPV);
0341         if (Eta < 0.)
0342           MPVsTECM->Fill(FitMPV);
0343         if (Eta > 0.)
0344           MPVsTECP->Fill(FitMPV);
0345         if (Thickness < 0.04) {
0346           MPV_Vs_EtaTECthin->Fill(Eta, FitMPV);
0347           MPV_Vs_PhiTECthin->Fill(Phi, FitMPV);
0348           MPVsTECthin->Fill(FitMPV);
0349           if (Eta > 0.)
0350             MPVsTECP1->Fill(FitMPV);
0351           if (Eta < 0.)
0352             MPVsTECM1->Fill(FitMPV);
0353         }
0354         if (Thickness > 0.04) {
0355           MPV_Vs_EtaTECthick->Fill(Eta, FitMPV);
0356           MPV_Vs_PhiTECthick->Fill(Phi, FitMPV);
0357           MPVsTECthick->Fill(FitMPV);
0358           if (Eta > 0.)
0359             MPVsTECP2->Fill(FitMPV);
0360           if (Eta < 0.)
0361             MPVsTECM2->Fill(FitMPV);
0362         }
0363       }
0364     }
0365 
0366     if (SubDet == 3 && PreviousGain != 0.)
0367       DiffWRTPrevGainTIB->Fill(Gain / PreviousGain);
0368     else if (SubDet == 4 && PreviousGain != 0.)
0369       DiffWRTPrevGainTID->Fill(Gain / PreviousGain);
0370     else if (SubDet == 5 && PreviousGain != 0.)
0371       DiffWRTPrevGainTOB->Fill(Gain / PreviousGain);
0372     else if (SubDet == 6 && PreviousGain != 0.)
0373       DiffWRTPrevGainTEC->Fill(Gain / PreviousGain);
0374 
0375     if (SubDet == 3)
0376       GainVsPrevGainTIB->Fill(PreviousGain, Gain);
0377     else if (SubDet == 4)
0378       GainVsPrevGainTID->Fill(PreviousGain, Gain);
0379     else if (SubDet == 5)
0380       GainVsPrevGainTOB->Fill(PreviousGain, Gain);
0381     else if (SubDet == 6)
0382       GainVsPrevGainTEC->Fill(PreviousGain, Gain);
0383   }
0384 }
0385 
0386 namespace {
0387 
0388   std::pair<double, double> findFitRange(TH1* inputHisto) {
0389     const auto prevErrorIgnoreLevel = gErrorIgnoreLevel;
0390     gErrorIgnoreLevel = kError;
0391     auto charge_clone = std::unique_ptr<TH1D>(dynamic_cast<TH1D*>(inputHisto->Rebin(10, "charge_clone")));
0392     gErrorIgnoreLevel = prevErrorIgnoreLevel;
0393     float max_content = -1;
0394     float xMax = -1;
0395     for (int i = 1; i < charge_clone->GetNbinsX() + 1; ++i) {
0396       const auto bin_content = charge_clone->GetBinContent(i);
0397       const auto bin_center = charge_clone->GetXaxis()->GetBinCenter(i);
0398       if ((bin_content > max_content) && (bin_center > 100.)) {
0399         max_content = bin_content;
0400         xMax = bin_center;
0401       }
0402     }
0403     return std::pair<double, double>(xMax - 100., xMax + 500.);
0404   }
0405 
0406   Double_t langaufun(Double_t* x, Double_t* par) {
0407     // Numeric constants
0408     Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0409     Double_t mpshift = -0.22278298;       // Landau maximum location
0410 
0411     // Control constants
0412     Double_t np = 100.0;  // number of convolution steps
0413     Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0414 
0415     // Variables
0416     Double_t xx;
0417     Double_t mpc;
0418     Double_t fland;
0419     Double_t sum = 0.0;
0420     Double_t xlow, xupp;
0421     Double_t step;
0422     Double_t i;
0423 
0424     // MP shift correction
0425     mpc = par[1] - mpshift * par[0];
0426 
0427     // Range of convolution integral
0428     xlow = x[0] - sc * par[3];
0429     xupp = x[0] + sc * par[3];
0430 
0431     step = (xupp - xlow) / np;
0432 
0433     // Convolution integral of Landau and Gaussian by sum
0434     for (i = 1.0; i <= np / 2; i++) {
0435       xx = xlow + (i - .5) * step;
0436       fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0437       sum += fland * TMath::Gaus(x[0], xx, par[3]);
0438 
0439       xx = xupp - (i - .5) * step;
0440       fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0441 
0442       sum += fland * TMath::Gaus(x[0], xx, par[3]);
0443     }
0444 
0445     return (par[2] * step * sum * invsq2pi / par[3]);
0446   }
0447 
0448   std::unique_ptr<TF1> langaufit(TH1D* his,
0449                                  Double_t* fitrange,
0450                                  Double_t* startvalues,
0451                                  Double_t* parlimitslo,
0452                                  Double_t* parlimitshi,
0453                                  Double_t* fitparams,
0454                                  Double_t* fiterrors,
0455                                  Double_t* ChiSqr,
0456                                  Int_t* NDF) {
0457     Int_t i;
0458     Char_t FunName[100];
0459 
0460     sprintf(FunName, "Fitfcn_%s", his->GetName());
0461 
0462     TF1* ffitold = dynamic_cast<TF1*>(gROOT->GetListOfFunctions()->FindObject(FunName));
0463     if (ffitold)
0464       delete ffitold;
0465 
0466     auto ffit = std::make_unique<TF1>(FunName, langaufun, fitrange[0], fitrange[1], 4);
0467     //TF1 *ffit = new TF1(FunName,langaufun,his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax(),4);
0468     ffit->SetParameters(startvalues);
0469     ffit->SetParNames("Width", "MP", "Area", "GSigma");
0470 
0471     for (i = 0; i < 4; i++) {
0472       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0473     }
0474 
0475     his->Fit(FunName, "QRB0");  // fit within specified range, use ParLimits, do not plot
0476 
0477     ffit->GetParameters(fitparams);  // obtain fit parameters
0478     for (i = 0; i < 4; i++) {
0479       fiterrors[i] = ffit->GetParError(i);  // obtain fit parameter errors
0480     }
0481     ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
0482     NDF[0] = ffit->GetNDF();           // obtain ndf
0483 
0484     return ffit;  // return fit function
0485   }
0486 }  // namespace
0487 
0488 //********************************************************************************//
0489 void SiStripGainsPCLHarvester::algoComputeMPVandGain(const MonitorElement* Charge_Vs_Index) {
0490   unsigned int I = 0;
0491   TH1D* Proj = nullptr;
0492   static constexpr double DEF_F = -9999.;
0493   double FitResults[6] = {DEF_F, DEF_F, DEF_F, DEF_F, DEF_F, DEF_F};
0494   double MPVmean = 300;
0495 
0496   if (Charge_Vs_Index == nullptr) {
0497     edm::LogError("SiStripGainsPCLHarvester")
0498         << "Harvesting: could not execute algoComputeMPVandGain method because " << m_calibrationMode
0499         << " statistics cannot be retrieved.\n"
0500         << "Please check if input contains " << m_calibrationMode << " data." << std::endl;
0501     return;
0502   }
0503 
0504   TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
0505 
0506   printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0507   printf("Fitting Charge Distribution  :");
0508   int TreeStep = APVsColl.size() / 50;
0509 
0510   for (auto it = APVsColl.begin(); it != APVsColl.end(); it++, I++) {
0511     if (I % TreeStep == 0) {
0512       printf(".");
0513       fflush(stdout);
0514     }
0515     std::shared_ptr<stAPVGain> APV = it->second;
0516     if (APV->Bin < 0)
0517       APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
0518 
0519     if (APV->isMasked) {
0520       APV->Gain = APV->PreviousGain;
0521       MASKED++;
0522       continue;
0523     }
0524 
0525     Proj = chvsidx->ProjectionY(
0526         "", chvsidx->GetXaxis()->FindBin(APV->Index), chvsidx->GetXaxis()->FindBin(APV->Index), "e");
0527     if (!Proj)
0528       continue;
0529 
0530     if (CalibrationLevel == 0) {
0531     } else if (CalibrationLevel == 1) {
0532       int SecondAPVId = APV->APVId;
0533       if (SecondAPVId % 2 == 0) {
0534         SecondAPVId = SecondAPVId + 1;
0535       } else {
0536         SecondAPVId = SecondAPVId - 1;
0537       }
0538       std::shared_ptr<stAPVGain> APV2 = APVsColl[(APV->DetId << 4) | SecondAPVId];
0539       if (APV2->Bin < 0)
0540         APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
0541       TH1D* Proj2 = chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e");
0542       if (Proj2) {
0543         Proj->Add(Proj2, 1);
0544         delete Proj2;
0545       }
0546     } else if (CalibrationLevel == 2) {
0547       for (unsigned int i = 0; i < 16; i++) {  //loop up to 6APV for Strip and up to 16 for Pixels
0548         auto tmpit = APVsColl.find((APV->DetId << 4) | i);
0549         if (tmpit == APVsColl.end())
0550           continue;
0551         std::shared_ptr<stAPVGain> APV2 = tmpit->second;
0552         if (APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)
0553           continue;
0554         if (APV2->Bin < 0)
0555           APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
0556         TH1D* Proj2 = chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e");
0557         if (Proj2) {
0558           Proj->Add(Proj2, 1);
0559           delete Proj2;
0560         }
0561       }
0562     } else {
0563       CalibrationLevel = 0;
0564       printf("Unknown Calibration Level, will assume %i\n", CalibrationLevel);
0565     }
0566 
0567     std::pair<double, double> fitRange{50., 5400.};
0568     if (fit_dataDrivenRange_) {
0569       fitRange = findFitRange(Proj);
0570     }
0571 
0572     const bool isTOBL5L6 =
0573         (DetId{APV->DetId}.subdetId() == StripSubdetector::TOB) && (tTopo_->tobLayer(APV->DetId) > 4);
0574     getPeakOfLandau(Proj,
0575                     FitResults,
0576                     fitRange.first,
0577                     fitRange.second,
0578                     (isTOBL5L6 ? fit_gaussianConvolutionTOBL56_ : fit_gaussianConvolution_));
0579 
0580     // throw if the fit results are not set
0581     assert(FitResults[0] != DEF_F);
0582 
0583     APV->FitMPV = FitResults[0];
0584     APV->FitMPVErr = FitResults[1];
0585     APV->FitWidth = FitResults[2];
0586     APV->FitWidthErr = FitResults[3];
0587     APV->FitChi2 = FitResults[4];
0588     APV->FitNorm = FitResults[5];
0589     APV->NEntries = Proj->GetEntries();
0590 
0591     // fall back to legacy fit in case of  very low chi2 probability
0592     if (APV->FitChi2 <= 0.1) {
0593       edm::LogInfo("SiStripGainsPCLHarvester")
0594           << "fit failed on this APV:" << APV->DetId << ":" << APV->APVId << " !" << std::endl;
0595 
0596       std::fill(FitResults, FitResults + 6, 0);
0597       fitRange = std::make_pair(50., 5400.);
0598 
0599       APV->FitGrade = fitgrade::B;
0600 
0601       getPeakOfLandau(Proj, FitResults, fitRange.first, fitRange.second, false);
0602 
0603       APV->FitMPV = FitResults[0];
0604       APV->FitMPVErr = FitResults[1];
0605       APV->FitWidth = FitResults[2];
0606       APV->FitWidthErr = FitResults[3];
0607       APV->FitChi2 = FitResults[4];
0608       APV->FitNorm = FitResults[5];
0609       APV->NEntries = Proj->GetEntries();
0610     } else {
0611       APV->FitGrade = fitgrade::A;
0612     }
0613 
0614     if (IsGoodLandauFit(FitResults)) {
0615       APV->Gain = APV->FitMPV / MPVmean;
0616       if (APV->SubDet > 2)
0617         GOOD++;
0618     } else {
0619       APV->Gain = APV->PreviousGain;
0620       if (APV->SubDet > 2)
0621         BAD++;
0622     }
0623     if (APV->Gain <= 0)
0624       APV->Gain = 1;
0625 
0626     delete Proj;
0627   }
0628   printf("\n");
0629 }
0630 
0631 //********************************************************************************//
0632 void SiStripGainsPCLHarvester::getPeakOfLandau(
0633     TH1* InputHisto, double* FitResults, double LowRange, double HighRange, bool gaussianConvolution) {
0634   // undo defaults (checked for assertion)
0635   FitResults[0] = -0.5;  //MPV
0636   FitResults[1] = 0;     //MPV error
0637   FitResults[2] = -0.5;  //Width
0638   FitResults[3] = 0;     //Width error
0639   FitResults[4] = -0.5;  //Fit Chi2/NDF
0640   FitResults[5] = 0;     //Normalization
0641 
0642   if (InputHisto->GetEntries() < MinNrEntries)
0643     return;
0644 
0645   if (gaussianConvolution) {
0646     // perform fit with landau convoluted with a gaussian
0647     Double_t fr[2] = {LowRange, HighRange};
0648     Double_t sv[4] = {25., 300., InputHisto->Integral(), 40.};
0649     Double_t pllo[4] = {0.5, 100., 1.0, 0.4};
0650     Double_t plhi[4] = {100., 500., 1.e7, 100.};
0651     Double_t fp[4], fpe[4];
0652     Double_t chisqr;
0653     Int_t ndf;
0654     auto fitsnr = langaufit(dynamic_cast<TH1D*>(InputHisto), fr, sv, pllo, plhi, fp, fpe, &chisqr, &ndf);
0655     FitResults[0] = fitsnr->GetMaximumX();  //MPV
0656     FitResults[1] = fpe[1];                 //MPV error // FIXME add error propagation
0657     FitResults[2] = fp[0];                  //Width
0658     FitResults[3] = fpe[0];                 //Width error
0659     FitResults[4] = chisqr / ndf;           //Fit Chi2/NDF
0660     FitResults[5] = fp[2];
0661   } else {
0662     // perform fit with standard landau
0663     TF1 MyLandau("MyLandau", "landau", LowRange, HighRange);
0664     MyLandau.SetParameter(1, 300);
0665     InputHisto->Fit(&MyLandau, "0QR WW");
0666     // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0667     FitResults[0] = MyLandau.GetParameter(1);                     //MPV
0668     FitResults[1] = MyLandau.GetParError(1);                      //MPV error
0669     FitResults[2] = MyLandau.GetParameter(2);                     //Width
0670     FitResults[3] = MyLandau.GetParError(2);                      //Width error
0671     FitResults[4] = MyLandau.GetChisquare() / MyLandau.GetNDF();  //Fit Chi2/NDF
0672     FitResults[5] = MyLandau.GetParameter(0);
0673   }
0674 }
0675 
0676 //********************************************************************************//
0677 bool SiStripGainsPCLHarvester::IsGoodLandauFit(double* FitResults) {
0678   if (FitResults[0] <= 0)
0679     return false;
0680   //   if(FitResults[1] > MaxMPVError   )return false;
0681   //   if(FitResults[4] > MaxChi2OverNDF)return false;
0682   return true;
0683 }
0684 
0685 //********************************************************************************//
0686 // ------------ method called once each job just before starting event loop  ------------
0687 void SiStripGainsPCLHarvester::checkBookAPVColls(const edm::EventSetup& es) {
0688   auto newBareTkGeomPtr = &es.getData(tkGeomToken_);
0689   if (newBareTkGeomPtr == bareTkGeomPtr_)
0690     return;  // already filled APVColls, nothing changed
0691 
0692   if (!bareTkGeomPtr_) {  // pointer not yet set: called the first time => fill the APVColls
0693     auto const& Det = newBareTkGeomPtr->dets();
0694 
0695     unsigned int Index = 0;
0696 
0697     for (unsigned int i = 0; i < Det.size(); i++) {
0698       DetId Detid = Det[i]->geographicalId();
0699       int SubDet = Detid.subdetId();
0700 
0701       if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
0702           SubDet == StripSubdetector::TEC) {
0703         auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
0704         if (!DetUnit)
0705           continue;
0706 
0707         const StripTopology& Topo = DetUnit->specificTopology();
0708         unsigned int NAPV = Topo.nstrips() / 128;
0709 
0710         for (unsigned int j = 0; j < NAPV; j++) {
0711           auto APV = std::make_shared<stAPVGain>();
0712           APV->Index = Index;
0713           APV->Bin = -1;
0714           APV->DetId = Detid.rawId();
0715           APV->APVId = j;
0716           APV->SubDet = SubDet;
0717           APV->FitMPV = -1;
0718           APV->FitMPVErr = -1;
0719           APV->FitWidth = -1;
0720           APV->FitWidthErr = -1;
0721           APV->FitChi2 = -1;
0722           APV->FitNorm = -1;
0723           APV->FitGrade = fitgrade::NONE;
0724           APV->Gain = -1;
0725           APV->PreviousGain = 1;
0726           APV->PreviousGainTick = 1;
0727           APV->x = DetUnit->position().basicVector().x();
0728           APV->y = DetUnit->position().basicVector().y();
0729           APV->z = DetUnit->position().basicVector().z();
0730           APV->Eta = DetUnit->position().basicVector().eta();
0731           APV->Phi = DetUnit->position().basicVector().phi();
0732           APV->R = DetUnit->position().basicVector().transverse();
0733           APV->Thickness = DetUnit->surface().bounds().thickness();
0734           APV->NEntries = 0;
0735           APV->isMasked = false;
0736 
0737           APVsCollOrdered.push_back(APV);
0738           APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0739           Index++;
0740           NStripAPVs++;
0741         }  // loop on APVs
0742       }    // if is Strips
0743     }      // loop on dets
0744 
0745     for (unsigned int i = 0; i < Det.size();
0746          i++) {  //Make two loop such that the Pixel information is added at the end --> make transition simpler
0747       DetId Detid = Det[i]->geographicalId();
0748       int SubDet = Detid.subdetId();
0749       if (SubDet == PixelSubdetector::PixelBarrel || SubDet == PixelSubdetector::PixelEndcap) {
0750         auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
0751         if (!DetUnit)
0752           continue;
0753 
0754         const PixelTopology& Topo = DetUnit->specificTopology();
0755         unsigned int NROCRow = Topo.nrows() / (80.);
0756         unsigned int NROCCol = Topo.ncolumns() / (52.);
0757 
0758         for (unsigned int j = 0; j < NROCRow; j++) {
0759           for (unsigned int i = 0; i < NROCCol; i++) {
0760             auto APV = std::make_shared<stAPVGain>();
0761             APV->Index = Index;
0762             APV->Bin = -1;
0763             APV->DetId = Detid.rawId();
0764             APV->APVId = (j << 3 | i);
0765             APV->SubDet = SubDet;
0766             APV->FitMPV = -1;
0767             APV->FitMPVErr = -1;
0768             APV->FitWidth = -1;
0769             APV->FitWidthErr = -1;
0770             APV->FitChi2 = -1;
0771             APV->FitGrade = fitgrade::NONE;
0772             APV->Gain = -1;
0773             APV->PreviousGain = 1;
0774             APV->PreviousGainTick = 1;
0775             APV->x = DetUnit->position().basicVector().x();
0776             APV->y = DetUnit->position().basicVector().y();
0777             APV->z = DetUnit->position().basicVector().z();
0778             APV->Eta = DetUnit->position().basicVector().eta();
0779             APV->Phi = DetUnit->position().basicVector().phi();
0780             APV->R = DetUnit->position().basicVector().transverse();
0781             APV->Thickness = DetUnit->surface().bounds().thickness();
0782             APV->isMasked = false;  //SiPixelQuality_->IsModuleBad(Detid.rawId());
0783             APV->NEntries = 0;
0784 
0785             APVsCollOrdered.push_back(APV);
0786             APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0787             Index++;
0788             NPixelDets++;
0789 
0790           }  // loop on ROC cols
0791         }    // loop on ROC rows
0792       }      // if Pixel
0793     }        // loop on Dets
0794   }          //if (!bareTkGeomPtr_) ...
0795   bareTkGeomPtr_ = newBareTkGeomPtr;
0796 }
0797 
0798 //********************************************************************************//
0799 bool SiStripGainsPCLHarvester::produceTagFilter(const MonitorElement* Charge_Vs_Index) {
0800   // The goal of this function is to check wether or not there is enough statistics
0801   // to produce a meaningful tag for the DB
0802 
0803   if (Charge_Vs_Index == nullptr) {
0804     edm::LogError("SiStripGainsPCLHarvester")
0805         << "produceTagFilter -> Return false: could not retrieve the " << m_calibrationMode << " statistics.\n"
0806         << "Please check if input contains " << m_calibrationMode << " data." << std::endl;
0807     return false;
0808   }
0809 
0810   float integral = (Charge_Vs_Index)->getTH2S()->Integral();
0811   if ((Charge_Vs_Index)->getTH2S()->Integral(0, NStripAPVs + 1, 0, 99999) < tagCondition_NClusters) {
0812     edm::LogWarning("SiStripGainsPCLHarvester")
0813         << "calibrationMode  -> " << m_calibrationMode << "\n"
0814         << "produceTagFilter -> Return false: Statistics is too low : " << integral << std::endl;
0815     return false;
0816   }
0817   if ((1.0 * GOOD) / (GOOD + BAD) < tagCondition_GoodFrac) {
0818     edm::LogWarning("SiStripGainsPCLHarvester")
0819         << "calibrationMode  -> " << m_calibrationMode << "\n"
0820         << "produceTagFilter ->  Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD + BAD)
0821         << std::endl;
0822     return false;
0823   }
0824   return true;
0825 }
0826 
0827 //********************************************************************************//
0828 std::unique_ptr<SiStripApvGain> SiStripGainsPCLHarvester::getNewObject(const MonitorElement* Charge_Vs_Index) {
0829   std::unique_ptr<SiStripApvGain> obj = std::make_unique<SiStripApvGain>();
0830 
0831   if (!produceTagFilter(Charge_Vs_Index)) {
0832     edm::LogWarning("SiStripGainsPCLHarvester")
0833         << "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << std::endl;
0834     return obj;
0835   } else {
0836     doStoreOnDB = true;
0837   }
0838 
0839   std::vector<float> theSiStripVector;
0840   unsigned int PreviousDetId = 0;
0841   for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
0842     std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
0843     if (APV == nullptr) {
0844       printf("Bug\n");
0845       continue;
0846     }
0847     if (APV->SubDet <= 2)
0848       continue;
0849     if (APV->DetId != PreviousDetId) {
0850       if (!theSiStripVector.empty()) {
0851         SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0852         if (!obj->put(PreviousDetId, range))
0853           printf("Bug to put detId = %i\n", PreviousDetId);
0854       }
0855       theSiStripVector.clear();
0856       PreviousDetId = APV->DetId;
0857     }
0858     theSiStripVector.push_back(APV->Gain);
0859 
0860     LogDebug("SiStripGainsPCLHarvester") << " DetId: " << APV->DetId << " APV:   " << APV->APVId
0861                                          << " Gain:  " << APV->Gain << std::endl;
0862   }
0863   if (!theSiStripVector.empty()) {
0864     SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0865     if (!obj->put(PreviousDetId, range))
0866       printf("Bug to put detId = %i\n", PreviousDetId);
0867   }
0868 
0869   return obj;
0870 }
0871 
0872 void SiStripGainsPCLHarvester::storeGainsTree(const TAxis* chVsIdxXaxis) const {
0873   unsigned int t_Index, t_Bin, t_DetId;
0874   unsigned char t_APVId, t_SubDet;
0875   float t_x, t_y, t_z, t_Eta, t_R, t_Phi, t_Thickness;
0876   float t_FitMPV, t_FitMPVErr, t_FitWidth, t_FitWidthErr, t_FitChi2NDF, t_FitNorm, t_FitGrade;
0877   double t_Gain, t_PrevGain, t_PrevGainTick, t_NEntries;
0878   bool t_isMasked;
0879   auto tree = edm::Service<TFileService>()->make<TTree>("APVGain", "APVGain");
0880   tree->Branch("Index", &t_Index, "Index/i");
0881   tree->Branch("Bin", &t_Bin, "Bin/i");
0882   tree->Branch("DetId", &t_DetId, "DetId/i");
0883   tree->Branch("APVId", &t_APVId, "APVId/b");
0884   tree->Branch("SubDet", &t_SubDet, "SubDet/b");
0885   tree->Branch("x", &t_x, "x/F");
0886   tree->Branch("y", &t_y, "y/F");
0887   tree->Branch("z", &t_z, "z/F");
0888   tree->Branch("Eta", &t_Eta, "Eta/F");
0889   tree->Branch("R", &t_R, "R/F");
0890   tree->Branch("Phi", &t_Phi, "Phi/F");
0891   tree->Branch("Thickness", &t_Thickness, "Thickness/F");
0892   tree->Branch("FitMPV", &t_FitMPV, "FitMPV/F");
0893   tree->Branch("FitMPVErr", &t_FitMPVErr, "FitMPVErr/F");
0894   tree->Branch("FitWidth", &t_FitWidth, "FitWidth/F");
0895   tree->Branch("FitWidthErr", &t_FitWidthErr, "FitWidthErr/F");
0896   tree->Branch("FitChi2NDF", &t_FitChi2NDF, "FitChi2NDF/F");
0897   tree->Branch("FitNorm", &t_FitNorm, "FitNorm/F");
0898   tree->Branch("FitGrade", &t_FitGrade, "FitGrade/F");
0899   tree->Branch("Gain", &t_Gain, "Gain/D");
0900   tree->Branch("PrevGain", &t_PrevGain, "PrevGain/D");
0901   tree->Branch("PrevGainTick", &t_PrevGainTick, "PrevGainTick/D");
0902   tree->Branch("NEntries", &t_NEntries, "NEntries/D");
0903   tree->Branch("isMasked", &t_isMasked, "isMasked/O");
0904 
0905   for (const auto& iAPV : APVsCollOrdered) {
0906     if (iAPV) {
0907       t_Index = iAPV->Index;
0908       t_Bin = chVsIdxXaxis->FindBin(iAPV->Index);
0909       t_DetId = iAPV->DetId;
0910       t_APVId = iAPV->APVId;
0911       t_SubDet = iAPV->SubDet;
0912       t_x = iAPV->x;
0913       t_y = iAPV->y;
0914       t_z = iAPV->z;
0915       t_Eta = iAPV->Eta;
0916       t_R = iAPV->R;
0917       t_Phi = iAPV->Phi;
0918       t_Thickness = iAPV->Thickness;
0919       t_FitMPV = iAPV->FitMPV;
0920       t_FitMPVErr = iAPV->FitMPVErr;
0921       t_FitWidth = iAPV->FitWidth;
0922       t_FitWidthErr = iAPV->FitWidthErr;
0923       t_FitChi2NDF = iAPV->FitChi2;
0924       t_FitNorm = iAPV->FitNorm;
0925       t_FitGrade = iAPV->FitGrade;
0926       t_Gain = iAPV->Gain;
0927       t_PrevGain = iAPV->PreviousGain;
0928       t_PrevGainTick = iAPV->PreviousGainTick;
0929       t_NEntries = iAPV->NEntries;
0930       t_isMasked = iAPV->isMasked;
0931 
0932       tree->Fill();
0933     }
0934   }
0935 }
0936 
0937 //********************************************************************************//
0938 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0939 void SiStripGainsPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0940   edm::ParameterSetDescription desc;
0941   desc.setUnknown();
0942   descriptions.addDefault(desc);
0943 }
0944 
0945 //********************************************************************************//
0946 void SiStripGainsPCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0947   if (!tTopo_) {
0948     tTopo_ = std::make_unique<TrackerTopology>(isetup.getData(tTopoToken_));
0949   }
0950 }