Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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