Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-06 07:38:12

0001 #include "CalibTracker/SiStripChannelGain/interface/SiStripGainsPCLWorker.h"
0002 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 
0005 #include <iostream>
0006 #include <sstream>
0007 
0008 //********************************************************************************//
0009 SiStripGainsPCLWorker::SiStripGainsPCLWorker(const edm::ParameterSet& iConfig) {
0010   MinTrackMomentum = iConfig.getUntrackedParameter<double>("minTrackMomentum", 3.0);
0011   MaxTrackMomentum = iConfig.getUntrackedParameter<double>("maxTrackMomentum", 99999.0);
0012   MinTrackEta = iConfig.getUntrackedParameter<double>("minTrackEta", -5.0);
0013   MaxTrackEta = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0);
0014   MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips", 2);
0015   MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits", 8);
0016   MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double>("MaxTrackChiOverNdf", 3);
0017   MaxTrackingIteration = iConfig.getUntrackedParameter<int>("MaxTrackingIteration", 7);
0018   AllowSaturation = iConfig.getUntrackedParameter<bool>("AllowSaturation", false);
0019   FirstSetOfConstants = iConfig.getUntrackedParameter<bool>("FirstSetOfConstants", true);
0020   Validation = iConfig.getUntrackedParameter<bool>("Validation", false);
0021   OldGainRemoving = iConfig.getUntrackedParameter<bool>("OldGainRemoving", false);
0022   useCalibration = iConfig.getUntrackedParameter<bool>("UseCalibration", false);
0023   doChargeMonitorPerPlane = iConfig.getUntrackedParameter<bool>("doChargeMonitorPerPlane", false);
0024   m_DQMdir = iConfig.getUntrackedParameter<std::string>("DQMdir", "AlCaReco/SiStripGains");
0025   m_calibrationMode = iConfig.getUntrackedParameter<std::string>("calibrationMode", "StdBunch");
0026   VChargeHisto = iConfig.getUntrackedParameter<std::vector<std::string>>("ChargeHisto");
0027 
0028   // fill in the mapping between the histogram indices and the (id,side,plane) tuple
0029   std::vector<std::pair<std::string, std::string>> hnames =
0030       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0031   for (unsigned int i = 0; i < hnames.size(); i++) {
0032     int id = APVGain::subdetectorId((hnames[i]).first);
0033     int side = APVGain::subdetectorSide((hnames[i]).first);
0034     int plane = APVGain::subdetectorPlane((hnames[i]).first);
0035     int thick = APVGain::thickness((hnames[i]).first);
0036     std::string s = hnames[i].first;
0037 
0038     auto loc = APVloc(thick, id, side, plane, s);
0039     theTopologyMap.insert(std::make_pair(i, loc));
0040   }
0041 
0042   //Set the monitoring element tag and store
0043   dqm_tag_.reserve(7);
0044   dqm_tag_.clear();
0045   dqm_tag_.push_back("StdBunch");    // statistic collection from Standard Collision Bunch @ 3.8 T
0046   dqm_tag_.push_back("StdBunch0T");  // statistic collection from Standard Collision Bunch @ 0 T
0047   dqm_tag_.push_back("AagBunch");    // statistic collection from First Collision After Abort Gap @ 3.8 T
0048   dqm_tag_.push_back("AagBunch0T");  // statistic collection from First Collision After Abort Gap @ 0 T
0049   dqm_tag_.push_back("IsoMuon");     // statistic collection from Isolated Muon @ 3.8 T
0050   dqm_tag_.push_back("IsoMuon0T");   // statistic collection from Isolated Muon @ 0 T
0051   dqm_tag_.push_back("Harvest");     // statistic collection: Harvest
0052 
0053   m_tracks_token = consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("tracks"));
0054   m_association_token = consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0055 
0056   tTopoToken_ = esConsumes();
0057   tkGeomTokenBR_ = esConsumes<edm::Transition::BeginRun>();
0058   tkGeomToken_ = esConsumes<>();
0059   gainToken_ = esConsumes<edm::Transition::BeginRun>();
0060   qualityToken_ = esConsumes<edm::Transition::BeginRun>();
0061 }
0062 
0063 //********************************************************************************//
0064 void SiStripGainsPCLWorker::dqmBeginRun(edm::Run const& run,
0065                                         edm::EventSetup const& iSetup,
0066                                         APVGain::APVGainHistograms& histograms) const {
0067   using namespace edm;
0068   static constexpr float defaultGainTick = 690. / 640.;
0069 
0070   // fills the APV collections at each begin run
0071   const TrackerGeometry* bareTkGeomPtr = &iSetup.getData(tkGeomTokenBR_);
0072   checkBookAPVColls(bareTkGeomPtr, histograms);
0073 
0074   const auto gainHandle = iSetup.getHandle(gainToken_);
0075   if (!gainHandle.isValid()) {
0076     edm::LogError("SiStripGainPCLWorker") << "gainHandle is not valid\n";
0077     exit(0);
0078   }
0079 
0080   const auto& siStripQuality = iSetup.getData(qualityToken_);
0081 
0082   for (unsigned int a = 0; a < histograms.APVsCollOrdered.size(); a++) {
0083     std::shared_ptr<stAPVGain> APV = histograms.APVsCollOrdered[a];
0084 
0085     if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
0086       continue;
0087 
0088     APV->isMasked = siStripQuality.IsApvBad(APV->DetId, APV->APVId);
0089 
0090     if (gainHandle->getNumberOfTags() != 2) {
0091       edm::LogError("SiStripGainPCLWorker") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
0092       fflush(stdout);
0093       exit(0);
0094     };
0095     float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
0096     if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
0097       edm::LogWarning("SiStripGainPCLWorker") << "WARNING: ParticleGain in the global tag changed\n";
0098     APV->PreviousGain = newPreviousGain;
0099 
0100     float newPreviousGainTick =
0101         APV->isMasked ? defaultGainTick : gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
0102     if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
0103       edm::LogWarning("SiStripGainPCLWorker")
0104           << "WARNING: TickMarkGain in the global tag changed\n"
0105           << std::endl
0106           << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
0107           << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
0108           << std::endl;
0109     }
0110     APV->PreviousGainTick = newPreviousGainTick;
0111   }
0112 }
0113 
0114 namespace {
0115   struct HitCluster {
0116     uint32_t det;
0117     const SiStripCluster* strip;
0118     const SiPixelCluster* pixel;
0119     HitCluster(uint32_t detId, const SiStripCluster* strip, const SiPixelCluster* pixel)
0120         : det(detId), strip(strip), pixel(pixel) {}
0121   };
0122   std::vector<HitCluster> getClusters(const TrackingRecHit* hit) {
0123     const auto simple1d = dynamic_cast<const SiStripRecHit1D*>(hit);
0124     const auto simple = dynamic_cast<const SiStripRecHit2D*>(hit);
0125     const auto matched = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0126     const auto pixel = dynamic_cast<const SiPixelRecHit*>(hit);
0127     std::vector<HitCluster> clusters;
0128     if (matched) {
0129       clusters.emplace_back(matched->monoId(), &matched->monoCluster(), nullptr);
0130       clusters.emplace_back(matched->stereoId(), &matched->stereoCluster(), nullptr);
0131     } else if (simple) {
0132       clusters.emplace_back(simple->geographicalId().rawId(), simple->cluster().get(), nullptr);
0133     } else if (simple1d) {
0134       clusters.emplace_back(simple1d->geographicalId().rawId(), simple1d->cluster().get(), nullptr);
0135     } else if (pixel) {
0136       clusters.emplace_back(pixel->geographicalId().rawId(), nullptr, pixel->cluster().get());
0137     }
0138     return clusters;
0139   }
0140 
0141   bool isFarFromBorder(const TrajectoryStateOnSurface& trajState, uint32_t detId, const TrackerGeometry* tGeom) {
0142     const auto gdu = tGeom->idToDetUnit(detId);
0143     if ((!dynamic_cast<const StripGeomDetUnit*>(gdu)) && (!dynamic_cast<const PixelGeomDetUnit*>(gdu))) {
0144       edm::LogWarning("SiStripGainCalibTableProducer")
0145           << "DetId " << detId << " does not seem to belong to the tracker";
0146       return false;
0147     }
0148     const auto plane = gdu->surface();
0149     const auto trapBounds = dynamic_cast<const TrapezoidalPlaneBounds*>(&plane.bounds());
0150     const auto rectBounds = dynamic_cast<const RectangularPlaneBounds*>(&plane.bounds());
0151 
0152     static constexpr double distFromBorder = 1.0;
0153     double halfLength = 0.;
0154     if (trapBounds) {
0155       halfLength = trapBounds->parameters()[3];
0156     } else if (rectBounds) {
0157       halfLength = .5 * gdu->surface().bounds().length();
0158     } else {
0159       return false;
0160     }
0161 
0162     const auto pos = trajState.localPosition();
0163     const auto posError = trajState.localError().positionError();
0164     if (std::abs(pos.y()) + posError.yy() >= (halfLength - distFromBorder))
0165       return false;
0166 
0167     return true;
0168   }
0169 }  // namespace
0170 
0171 //********************************************************************************//
0172 // ------------ method called for each event  ------------
0173 void SiStripGainsPCLWorker::dqmAnalyze(edm::Event const& iEvent,
0174                                        edm::EventSetup const& iSetup,
0175                                        APVGain::APVGainHistograms const& histograms) const {
0176   using namespace edm;
0177 
0178   unsigned int eventnumber = iEvent.id().event();
0179   unsigned int runnumber = iEvent.id().run();
0180 
0181   edm::LogInfo("SiStripGainsPCLWorker") << "Processing run " << runnumber << " and event " << eventnumber << std::endl;
0182 
0183   const TrackerTopology* topo = &iSetup.getData(tTopoToken_);
0184   const TrackerGeometry* tGeom = &iSetup.getData(tkGeomToken_);
0185 
0186   // Event data handles
0187   edm::Handle<edm::View<reco::Track>> tracks;
0188   iEvent.getByToken(m_tracks_token, tracks);
0189   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociations;
0190   iEvent.getByToken(m_association_token, trajTrackAssociations);
0191 
0192   for (const auto& elem : theTopologyMap) {
0193     LogDebug("SiStripGainsPCLWorker") << elem.first << " - " << elem.second.m_string << " "
0194                                       << elem.second.m_subdetectorId << " " << elem.second.m_subdetectorSide << " "
0195                                       << elem.second.m_subdetectorPlane << std::endl;
0196   }
0197 
0198   LogDebug("SiStripGainsPCLWorker") << "for mode" << m_calibrationMode << std::endl;
0199 
0200   int elepos = statCollectionFromMode(m_calibrationMode.c_str());
0201 
0202   std::size_t nStoredClusters{0};
0203   for (const auto& assoc : *trajTrackAssociations) {
0204     const auto traj = assoc.key.get();
0205     const auto track = assoc.val.get();
0206 
0207     if ((track->eta() < MinTrackEta) || (track->eta() > MaxTrackEta) || (track->p() < MinTrackMomentum) ||
0208         (track->p() > MaxTrackMomentum) || (track->numberOfValidHits() < MinTrackHits) ||
0209         ((track->chi2() / track->ndof()) > MaxTrackChiOverNdf) || (track->algo() > MaxTrackingIteration))
0210       continue;
0211 
0212     int iCluster{-1};
0213     for (const auto& meas : traj->measurements()) {
0214       const auto& trajState = meas.updatedState();
0215       if (!trajState.isValid())
0216         continue;
0217 
0218       // there can be 2 (stereo module), 1 (no stereo module), or 0 (no pixel or strip hit) clusters
0219       auto clusters = getClusters(meas.recHit()->hit());
0220       for (const auto hitCluster : clusters) {
0221         ++iCluster;
0222         bool saturation = false;
0223         bool overlapping = false;
0224         unsigned int charge = 0;
0225         int firstStrip = 0;
0226         unsigned int nStrips = 0;
0227         if (hitCluster.strip) {
0228           const auto& ampls = hitCluster.strip->amplitudes();
0229           firstStrip = hitCluster.strip->firstStrip();
0230           nStrips = ampls.size();
0231           charge = hitCluster.strip->charge();
0232           saturation = std::any_of(ampls.begin(), ampls.end(), [](uint8_t amp) { return amp >= 254; });
0233 
0234           overlapping = (((firstStrip % 128) == 0) || ((firstStrip / 128) != ((firstStrip + int(nStrips)) / 128)));
0235         } else if (hitCluster.pixel) {
0236           const auto& ampls = hitCluster.pixel->pixelADC();
0237           const int firstRow = hitCluster.pixel->minPixelRow();
0238           const int firstCol = hitCluster.pixel->minPixelCol();
0239           firstStrip = ((firstRow / 80) << 3 | (firstCol / 52)) * 128;  //Hack to save the APVId
0240           nStrips = 0;
0241           for (const auto amp : ampls) {
0242             charge += amp;
0243             if (amp >= 254)
0244               saturation = true;
0245           }
0246         }
0247         // works for both strip and pixel thanks to firstStrip encoding for pixel above, as in the calibTree
0248         std::shared_ptr<stAPVGain> APV = histograms.APVsColl.at((hitCluster.det << 4) | (firstStrip / 128));
0249 
0250         const auto farFromEdge = (hitCluster.strip ? isFarFromBorder(trajState, hitCluster.det, tGeom) : true);
0251         if ((APV->SubDet > 2) &&
0252             ((!farFromEdge) || overlapping || (saturation && !AllowSaturation) || (nStrips > MaxNrStrips)))
0253           continue;
0254 
0255         int clusterCharge = 0;
0256         if (APV->SubDet > 2) {  // strip
0257           if (useCalibration || !FirstSetOfConstants) {
0258             saturation = false;
0259             for (const auto origCharge : hitCluster.strip->amplitudes()) {
0260               int stripCharge;
0261               if (useCalibration) {
0262                 if (FirstSetOfConstants) {
0263                   stripCharge = int(origCharge / APV->CalibGain);
0264                 } else {
0265                   stripCharge = int(origCharge * (APV->PreviousGain / APV->CalibGain));
0266                 }
0267               } else {
0268                 if (FirstSetOfConstants) {
0269                   stripCharge = origCharge;
0270                 } else {
0271                   stripCharge = int(origCharge * APV->PreviousGain);
0272                 }
0273               }
0274               if (stripCharge > 1024) {
0275                 stripCharge = 255;
0276                 saturation = true;
0277               } else if (stripCharge > 254) {
0278                 stripCharge = 254;
0279                 saturation = true;
0280               }
0281               clusterCharge += stripCharge;
0282             }
0283             if (saturation && !AllowSaturation)
0284               continue;
0285           } else {
0286             clusterCharge = charge;
0287           }
0288         } else {                           // pixel
0289           clusterCharge = charge / 265.0;  //expected scale factor between pixel and strip charge
0290         }
0291 
0292         const auto trackDir = trajState.localDirection();
0293         const auto path = (10. * APV->Thickness) / std::abs(trackDir.z() / trackDir.mag());
0294         double ClusterChargeOverPath = ((double)clusterCharge) / path;
0295         if (APV->SubDet > 2) {
0296           if (Validation) {
0297             ClusterChargeOverPath /= APV->PreviousGain;
0298           }
0299           if (OldGainRemoving) {
0300             ClusterChargeOverPath *= APV->PreviousGain;
0301           }
0302         } else {
0303           // keep processing of pixel cluster charge until here
0304           continue;
0305         }
0306         ++nStoredClusters;
0307 
0308         // real histogram for calibration
0309         histograms.Charge_Vs_Index[elepos]->Fill(APV->Index, ClusterChargeOverPath);
0310         LogDebug("SiStripGainsPCLWorker")
0311             << " for mode " << m_calibrationMode << "\n"
0312             << " i " << iCluster << " useCalibration " << useCalibration << " FirstSetOfConstants "
0313             << FirstSetOfConstants << " APV->PreviousGain " << APV->PreviousGain << " APV->CalibGain " << APV->CalibGain
0314             << " APV->DetId " << APV->DetId << " APV->Index " << APV->Index << " Charge " << clusterCharge << " Path "
0315             << path << " ClusterChargeOverPath " << ClusterChargeOverPath << std::endl;
0316 
0317         // Fill monitoring histograms
0318         int mCharge1 = 0;
0319         for (const auto sCharge : hitCluster.strip->amplitudes()) {
0320           if (sCharge > 254) {
0321             mCharge1 += 254;
0322           } else {
0323             mCharge1 += sCharge;
0324           }
0325         }
0326         // Revome gains for monitoring
0327         int mCharge2 = mCharge1 * APV->PreviousGain;                          // remove G2
0328         int mCharge3 = mCharge1 * APV->PreviousGainTick;                      // remove G1
0329         int mCharge4 = mCharge1 * APV->PreviousGain * APV->PreviousGainTick;  // remove G1 and G2
0330 
0331         LogDebug("SiStripGainsPCLWorker") << " full charge " << mCharge1 << " remove G2 " << mCharge2 << " remove G1 "
0332                                           << mCharge3 << " remove G1*G2 " << mCharge4 << std::endl;
0333 
0334         auto indices = APVGain::FetchIndices(theTopologyMap, hitCluster.det, topo);
0335 
0336         for (auto m : indices)
0337           histograms.Charge_1[elepos][m]->Fill(((double)mCharge1) / path);
0338         for (auto m : indices)
0339           histograms.Charge_2[elepos][m]->Fill(((double)mCharge2) / path);
0340         for (auto m : indices)
0341           histograms.Charge_3[elepos][m]->Fill(((double)mCharge3) / path);
0342         for (auto m : indices)
0343           histograms.Charge_4[elepos][m]->Fill(((double)mCharge4) / path);
0344 
0345         if (APV->SubDet == StripSubdetector::TIB) {
0346           histograms.Charge_Vs_PathlengthTIB[elepos]->Fill(path, clusterCharge);  // TIB
0347 
0348         } else if (APV->SubDet == StripSubdetector::TOB) {
0349           histograms.Charge_Vs_PathlengthTOB[elepos]->Fill(path, clusterCharge);  // TOB
0350 
0351         } else if (APV->SubDet == StripSubdetector::TID) {
0352           if (APV->Eta < 0) {
0353             histograms.Charge_Vs_PathlengthTIDM[elepos]->Fill(path, clusterCharge);
0354           }  // TID minus
0355           else if (APV->Eta > 0) {
0356             histograms.Charge_Vs_PathlengthTIDP[elepos]->Fill(path, clusterCharge);
0357           }  // TID plus
0358 
0359         } else if (APV->SubDet == StripSubdetector::TEC) {
0360           if (APV->Eta < 0) {
0361             if (APV->Thickness < 0.04) {
0362               histograms.Charge_Vs_PathlengthTECM1[elepos]->Fill(path, clusterCharge);
0363             }  // TEC minus, type 1
0364             else if (APV->Thickness > 0.04) {
0365               histograms.Charge_Vs_PathlengthTECM2[elepos]->Fill(path, clusterCharge);
0366             }  // TEC minus, type 2
0367           } else if (APV->Eta > 0) {
0368             if (APV->Thickness < 0.04) {
0369               histograms.Charge_Vs_PathlengthTECP1[elepos]->Fill(path, clusterCharge);
0370             }  // TEC plus, type 1
0371             else if (APV->Thickness > 0.04) {
0372               histograms.Charge_Vs_PathlengthTECP2[elepos]->Fill(path, clusterCharge);
0373             }  // TEC plus, type 2
0374           }
0375         }
0376       }
0377     }
0378   }
0379 
0380   histograms.EventStats->Fill(0., 0., 1);
0381   histograms.EventStats->Fill(1., 0., tracks->size());
0382   histograms.EventStats->Fill(2., 0., nStoredClusters);
0383 
0384   //LogDebug("SiStripGainsPCLWorker")<<" for mode"<< m_calibrationMode
0385   //                   <<" entries in histogram:"<< histograms.Charge_Vs_Index[elepos].getEntries()
0386   //                   <<std::endl;
0387 }
0388 
0389 //********************************************************************************//
0390 // ------------ method called once each job just before starting event loop  ------------
0391 void SiStripGainsPCLWorker::checkBookAPVColls(const TrackerGeometry* bareTkGeomPtr,
0392                                               APVGain::APVGainHistograms& histograms) const {
0393   if (bareTkGeomPtr) {  // pointer not yet set: called the first time => fill the APVColls
0394     auto const& Det = bareTkGeomPtr->dets();
0395 
0396     edm::LogInfo("SiStripGainsPCLWorker") << " Resetting APV struct" << std::endl;
0397 
0398     unsigned int Index = 0;
0399 
0400     for (unsigned int i = 0; i < Det.size(); i++) {
0401       DetId Detid = Det[i]->geographicalId();
0402       int SubDet = Detid.subdetId();
0403 
0404       if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
0405           SubDet == StripSubdetector::TEC) {
0406         auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
0407         if (!DetUnit)
0408           continue;
0409 
0410         const StripTopology& Topo = DetUnit->specificTopology();
0411         unsigned int NAPV = Topo.nstrips() / 128;
0412 
0413         for (unsigned int j = 0; j < NAPV; j++) {
0414           auto APV = std::make_shared<stAPVGain>();
0415           APV->Index = Index;
0416           APV->Bin = -1;
0417           APV->DetId = Detid.rawId();
0418           APV->APVId = j;
0419           APV->SubDet = SubDet;
0420           APV->FitMPV = -1;
0421           APV->FitMPVErr = -1;
0422           APV->FitWidth = -1;
0423           APV->FitWidthErr = -1;
0424           APV->FitChi2 = -1;
0425           APV->FitNorm = -1;
0426           APV->Gain = -1;
0427           APV->PreviousGain = 1;
0428           APV->PreviousGainTick = 1;
0429           APV->x = DetUnit->position().basicVector().x();
0430           APV->y = DetUnit->position().basicVector().y();
0431           APV->z = DetUnit->position().basicVector().z();
0432           APV->Eta = DetUnit->position().basicVector().eta();
0433           APV->Phi = DetUnit->position().basicVector().phi();
0434           APV->R = DetUnit->position().basicVector().transverse();
0435           APV->Thickness = DetUnit->surface().bounds().thickness();
0436           APV->NEntries = 0;
0437           APV->isMasked = false;
0438 
0439           histograms.APVsCollOrdered.push_back(APV);
0440           histograms.APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0441           Index++;
0442           histograms.NStripAPVs++;
0443         }  // loop on APVs
0444       }    // if is Strips
0445     }      // loop on dets
0446 
0447     for (unsigned int i = 0; i < Det.size();
0448          i++) {  //Make two loop such that the Pixel information is added at the end --> make transition simpler
0449       DetId Detid = Det[i]->geographicalId();
0450       int SubDet = Detid.subdetId();
0451       if (SubDet == PixelSubdetector::PixelBarrel || SubDet == PixelSubdetector::PixelEndcap) {
0452         auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
0453         if (!DetUnit)
0454           continue;
0455 
0456         const PixelTopology& Topo = DetUnit->specificTopology();
0457         unsigned int NROCRow = Topo.nrows() / (80.);
0458         unsigned int NROCCol = Topo.ncolumns() / (52.);
0459 
0460         for (unsigned int j = 0; j < NROCRow; j++) {
0461           for (unsigned int i = 0; i < NROCCol; i++) {
0462             auto APV = std::make_shared<stAPVGain>();
0463             APV->Index = Index;
0464             APV->Bin = -1;
0465             APV->DetId = Detid.rawId();
0466             APV->APVId = (j << 3 | i);
0467             APV->SubDet = SubDet;
0468             APV->FitMPV = -1;
0469             APV->FitMPVErr = -1;
0470             APV->FitWidth = -1;
0471             APV->FitWidthErr = -1;
0472             APV->FitChi2 = -1;
0473             APV->Gain = -1;
0474             APV->PreviousGain = 1;
0475             APV->PreviousGainTick = 1;
0476             APV->x = DetUnit->position().basicVector().x();
0477             APV->y = DetUnit->position().basicVector().y();
0478             APV->z = DetUnit->position().basicVector().z();
0479             APV->Eta = DetUnit->position().basicVector().eta();
0480             APV->Phi = DetUnit->position().basicVector().phi();
0481             APV->R = DetUnit->position().basicVector().transverse();
0482             APV->Thickness = DetUnit->surface().bounds().thickness();
0483             APV->isMasked = false;  //SiPixelQuality_->IsModuleBad(Detid.rawId());
0484             APV->NEntries = 0;
0485 
0486             histograms.APVsCollOrdered.push_back(APV);
0487             histograms.APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
0488             Index++;
0489             histograms.NPixelDets++;
0490 
0491           }  // loop on ROC cols
0492         }    // loop on ROC rows
0493       }      // if Pixel
0494     }        // loop on Dets
0495   }          //if (!bareTkGeomPtr_) ...
0496 }
0497 
0498 //********************************************************************************//
0499 void SiStripGainsPCLWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0500   edm::ParameterSetDescription desc;
0501   desc.setUnknown();
0502   descriptions.addDefault(desc);
0503 }
0504 
0505 //********************************************************************************//
0506 void SiStripGainsPCLWorker::bookHistograms(DQMStore::IBooker& ibooker,
0507                                            edm::Run const& run,
0508                                            edm::EventSetup const& setup,
0509                                            APVGain::APVGainHistograms& histograms) const {
0510   ibooker.cd();
0511   std::string dqm_dir = m_DQMdir;
0512   const char* tag = dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str();
0513 
0514   edm::LogInfo("SiStripGainsPCLWorker") << "Setting " << dqm_dir << " in DQM and booking histograms for tag " << tag
0515                                         << std::endl;
0516 
0517   ibooker.setCurrentFolder(dqm_dir);
0518 
0519   // this MonitorElement is created to log the number of events / tracks and clusters used
0520   // by the calibration algorithm
0521 
0522   histograms.EventStats = ibooker.book2I("EventStats", "Statistics", 3, -0.5, 2.5, 1, 0, 1);
0523   histograms.EventStats->setBinLabel(1, "events count", 1);
0524   histograms.EventStats->setBinLabel(2, "tracks count", 1);
0525   histograms.EventStats->setBinLabel(3, "clusters count", 1);
0526 
0527   std::string stag(tag);
0528   if (!stag.empty() && stag[0] != '_')
0529     stag.insert(0, 1, '_');
0530 
0531   std::string cvi = std::string("Charge_Vs_Index") + stag;
0532   std::string cvpTIB = std::string("Charge_Vs_PathlengthTIB") + stag;
0533   std::string cvpTOB = std::string("Charge_Vs_PathlengthTOB") + stag;
0534   std::string cvpTIDP = std::string("Charge_Vs_PathlengthTIDP") + stag;
0535   std::string cvpTIDM = std::string("Charge_Vs_PathlengthTIDM") + stag;
0536   std::string cvpTECP1 = std::string("Charge_Vs_PathlengthTECP1") + stag;
0537   std::string cvpTECP2 = std::string("Charge_Vs_PathlengthTECP2") + stag;
0538   std::string cvpTECM1 = std::string("Charge_Vs_PathlengthTECM1") + stag;
0539   std::string cvpTECM2 = std::string("Charge_Vs_PathlengthTECM2") + stag;
0540 
0541   int elepos = statCollectionFromMode(tag);
0542 
0543   histograms.Charge_Vs_Index.reserve(dqm_tag_.size());
0544   histograms.Charge_Vs_PathlengthTIB.reserve(dqm_tag_.size());
0545   histograms.Charge_Vs_PathlengthTOB.reserve(dqm_tag_.size());
0546   histograms.Charge_Vs_PathlengthTIDP.reserve(dqm_tag_.size());
0547   histograms.Charge_Vs_PathlengthTIDM.reserve(dqm_tag_.size());
0548   histograms.Charge_Vs_PathlengthTECP1.reserve(dqm_tag_.size());
0549   histograms.Charge_Vs_PathlengthTECP2.reserve(dqm_tag_.size());
0550   histograms.Charge_Vs_PathlengthTECM1.reserve(dqm_tag_.size());
0551   histograms.Charge_Vs_PathlengthTECM2.reserve(dqm_tag_.size());
0552 
0553   // The cluster charge is stored by exploiting a non uniform binning in order
0554   // reduce the histogram memory size. The bin width is relaxed with a falling
0555   // exponential function and the bin boundaries are stored in the binYarray.
0556   // The binXarray is used to provide as many bins as the APVs.
0557   //
0558   // More details about this implementations are here:
0559   // https://indico.cern.ch/event/649344/contributions/2672267/attachments/1498323/2332518/OptimizeChHisto.pdf
0560 
0561   std::vector<float> binXarray;
0562   binXarray.reserve(histograms.NStripAPVs + 1);
0563   for (unsigned int a = 0; a <= histograms.NStripAPVs; a++) {
0564     binXarray.push_back((float)a);
0565   }
0566 
0567   std::array<float, 688> binYarray;
0568   double p0 = 5.445;
0569   double p1 = 0.002113;
0570   double p2 = 69.01576;
0571   double y = 0.;
0572   for (int b = 0; b < 687; b++) {
0573     binYarray[b] = y;
0574     if (y <= 902.)
0575       y = y + 2.;
0576     else
0577       y = (p0 - log(exp(p0 - p1 * y) - p2 * p1)) / p1;
0578   }
0579   binYarray[687] = 4000.;
0580 
0581   histograms.Charge_1[elepos].clear();
0582   histograms.Charge_2[elepos].clear();
0583   histograms.Charge_3[elepos].clear();
0584   histograms.Charge_4[elepos].clear();
0585 
0586   auto it = histograms.Charge_Vs_Index.begin();
0587   histograms.Charge_Vs_Index.insert(
0588       it + elepos,
0589       ibooker.book2S(cvi.c_str(), cvi.c_str(), histograms.NStripAPVs, &binXarray[0], 687, binYarray.data()));
0590 
0591   it = histograms.Charge_Vs_PathlengthTIB.begin();
0592   histograms.Charge_Vs_PathlengthTIB.insert(it + elepos,
0593                                             ibooker.book2S(cvpTIB.c_str(), cvpTIB.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0594 
0595   it = histograms.Charge_Vs_PathlengthTOB.begin();
0596   histograms.Charge_Vs_PathlengthTOB.insert(it + elepos,
0597                                             ibooker.book2S(cvpTOB.c_str(), cvpTOB.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0598 
0599   it = histograms.Charge_Vs_PathlengthTIDP.begin();
0600   histograms.Charge_Vs_PathlengthTIDP.insert(
0601       it + elepos, ibooker.book2S(cvpTIDP.c_str(), cvpTIDP.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0602 
0603   it = histograms.Charge_Vs_PathlengthTIDM.begin();
0604   histograms.Charge_Vs_PathlengthTIDM.insert(
0605       it + elepos, ibooker.book2S(cvpTIDM.c_str(), cvpTIDM.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0606 
0607   it = histograms.Charge_Vs_PathlengthTECP1.begin();
0608   histograms.Charge_Vs_PathlengthTECP1.insert(
0609       it + elepos, ibooker.book2S(cvpTECP1.c_str(), cvpTECP1.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0610 
0611   it = histograms.Charge_Vs_PathlengthTECP2.begin();
0612   histograms.Charge_Vs_PathlengthTECP2.insert(
0613       it + elepos, ibooker.book2S(cvpTECP2.c_str(), cvpTECP2.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0614 
0615   it = histograms.Charge_Vs_PathlengthTECM1.begin();
0616   histograms.Charge_Vs_PathlengthTECM1.insert(
0617       it + elepos, ibooker.book2S(cvpTECM1.c_str(), cvpTECM1.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0618 
0619   it = histograms.Charge_Vs_PathlengthTECM2.begin();
0620   histograms.Charge_Vs_PathlengthTECM2.insert(
0621       it + elepos, ibooker.book2S(cvpTECM2.c_str(), cvpTECM2.c_str(), 20, 0.3, 1.3, 250, 0, 2000));
0622 
0623   std::vector<std::pair<std::string, std::string>> hnames =
0624       APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "");
0625   for (unsigned int i = 0; i < hnames.size(); i++) {
0626     std::string htag = (hnames[i]).first + stag;
0627     histograms.Charge_1[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0628   }
0629 
0630   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG2");
0631   for (unsigned int i = 0; i < hnames.size(); i++) {
0632     std::string htag = (hnames[i]).first + stag;
0633     histograms.Charge_2[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0634   }
0635 
0636   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1");
0637   for (unsigned int i = 0; i < hnames.size(); i++) {
0638     std::string htag = (hnames[i]).first + stag;
0639     histograms.Charge_3[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0640   }
0641 
0642   hnames = APVGain::monHnames(VChargeHisto, doChargeMonitorPerPlane, "woG1G2");
0643   for (unsigned int i = 0; i < hnames.size(); i++) {
0644     std::string htag = (hnames[i]).first + stag;
0645     histograms.Charge_4[elepos].push_back(ibooker.book1DD(htag.c_str(), (hnames[i]).second.c_str(), 100, 0., 1000.));
0646   }
0647 }