Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!