Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-22 02:24:11 

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0006 #include "FWCore/ParameterSet/interface/FileInPath.h"
0007 
0008 #include "DataFormats/TrackReco/interface/Track.h"
0009 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
0010 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0011 #include "DataFormats/TrackReco/interface/DeDxHitInfo.h"
0012 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0013 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h"
0014 #include "CondFormats/PhysicsToolsObjects/interface/DeDxCalibration.h"
0015 #include "CondFormats/DataRecord/interface/DeDxCalibrationRcd.h"
0016 #include "CondFormats/SiPixelObjects/interface/PixelIndices.h"
0017 
0018 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0019 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0020 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0023 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineService.h"
0024 
0025 #include <fstream>
0026 
0027 class DeDxHitCalibrator : public edm::stream::EDProducer<> {
0028 public:
0029   static constexpr int kIsNormal = 0, kIsBelow = 1, kIsOver = 2;
0030   static constexpr int PXB = 0, PXF = 1, TIB = 2, TID = 3, TOB = 4, TECThin = 5, TECThick = 6;
0031   typedef std::pair<uint32_t, unsigned char> ChipId;
0032 
0033   explicit DeDxHitCalibrator(const edm::ParameterSet&);
0034   ~DeDxHitCalibrator() override = default;
0035   static void fillDescriptions(edm::ConfigurationDescriptions&);
0036 
0037 private:
0038   void beginRun(edm::Run const&, const edm::EventSetup&) override;
0039   void produce(edm::Event&, const edm::EventSetup&) override;
0040 
0041   int getDetId(const DetId&, const float&);
0042   float correctEnergy(const float&, const ChipId&);
0043   void processHitInfo(const reco::DeDxHitInfo&,
0044                       const float& trackMomentum,
0045                       reco::DeDxHitCollection&,
0046                       reco::DeDxHitCollection&);
0047 
0048   double getChi2(const std::vector<double>&, const std::vector<std::pair<double, int> >&, const double&, const double&);
0049   void getAlphaBeta(const std::vector<double>&,
0050                     const std::vector<std::pair<double, int> >&,
0051                     CLHEP::HepMatrix&,
0052                     CLHEP::HepVector&,
0053                     const std::vector<bool>&,
0054                     const double&,
0055                     const double&);
0056   std::pair<double, double> fitStripCluster(const std::vector<std::pair<double, int> >&, const double&, const double&);
0057 
0058   const bool applyGain_;
0059   const double MeVPerElectron_;
0060   const int VCaltoElectronGain_, VCaltoElectronGain_L1_, VCaltoElectronOffset_, VCaltoElectronOffset_L1_;
0061   const int pixelSaturationThr_;
0062   const edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0063   const edm::EDGetTokenT<reco::DeDxHitInfoAss> dedxHitInfoToken_;
0064   const edm::ESGetToken<DeDxCalibration, DeDxCalibrationRcd> dedxCalibToken_;
0065   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0066   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tkTopoToken_;
0067 
0068   SiPixelGainCalibrationOfflineService pixelCalib_;
0069   edm::ESHandle<DeDxCalibration> dedxCalib_;
0070   edm::ESHandle<TrackerGeometry> tkGeom_;
0071   edm::ESHandle<TrackerTopology> tkTopo_;
0072 };
0073 
0074 DeDxHitCalibrator::DeDxHitCalibrator(const edm::ParameterSet& iConfig)
0075     : applyGain_(iConfig.getParameter<bool>("applyGain")),
0076       MeVPerElectron_(iConfig.getParameter<double>("MeVPerElectron")),
0077       VCaltoElectronGain_(iConfig.getParameter<int>("VCaltoElectronGain")),
0078       VCaltoElectronGain_L1_(iConfig.getParameter<int>("VCaltoElectronGain_L1")),
0079       VCaltoElectronOffset_(iConfig.getParameter<int>("VCaltoElectronOffset")),
0080       VCaltoElectronOffset_L1_(iConfig.getParameter<int>("VCaltoElectronOffset_L1")),
0081       pixelSaturationThr_(iConfig.getParameter<int>("pixelSaturationThr")),
0082       tracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trackProducer"))),
0083       dedxHitInfoToken_(consumes<reco::DeDxHitInfoAss>(iConfig.getParameter<edm::InputTag>("dedxHitInfo"))),
0084       dedxCalibToken_(esConsumes<DeDxCalibration, DeDxCalibrationRcd, edm::Transition::BeginRun>()),
0085       tkGeomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0086       tkTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0087       pixelCalib_(iConfig, consumesCollector()) {
0088   produces<reco::TrackDeDxHitsCollection>("PixelHits");
0089   produces<reco::TrackDeDxHitsCollection>("StripHits");
0090 }
0091 
0092 void DeDxHitCalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0093   edm::ParameterSetDescription desc;
0094   desc.add<bool>("applyGain", true);
0095   desc.add<double>("MeVPerElectron", 3.61e-06);
0096   desc.add<int>("VCaltoElectronGain", 65);
0097   desc.add<int>("VCaltoElectronGain_L1", 65);
0098   desc.add<int>("VCaltoElectronOffset", -414);
0099   desc.add<int>("VCaltoElectronOffset_L1", -414);
0100   desc.add<int>("pixelSaturationThr", 254);
0101   desc.add<edm::InputTag>("trackProducer", edm::InputTag("generalTracks"));
0102   desc.add<edm::InputTag>("dedxHitInfo", edm::InputTag("dedxHitInfo"));
0103   descriptions.add("dedxHitCalibrator", desc);
0104 }
0105 
0106 void DeDxHitCalibrator::beginRun(edm::Run const&, const edm::EventSetup& iSetup) {
0107   dedxCalib_ = iSetup.getHandle(dedxCalibToken_);
0108   tkGeom_ = iSetup.getHandle(tkGeomToken_);
0109   tkTopo_ = iSetup.getHandle(tkTopoToken_);
0110 }
0111 
0112 void DeDxHitCalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113   const auto& tracks = iEvent.getHandle(tracksToken_);
0114   const auto& dedxHitInfo = iEvent.get(dedxHitInfoToken_);
0115   pixelCalib_.setESObjects(iSetup);
0116 
0117   // creates the output collection
0118   auto pixelTrackDeDxHitAss = std::make_unique<reco::TrackDeDxHitsCollection>(reco::TrackRefProd(tracks));
0119   auto stripTrackDeDxHitAss = std::make_unique<reco::TrackDeDxHitsCollection>(reco::TrackRefProd(tracks));
0120 
0121   for (size_t i = 0; i < tracks->size(); i++) {
0122     const auto& track = reco::TrackRef(tracks, i);
0123     const auto& dedxHits = dedxHitInfo[track];
0124     if (dedxHits.isNull())
0125       continue;
0126     reco::DeDxHitCollection pixelHits, stripHits;
0127     processHitInfo(*dedxHits, track->p(), pixelHits, stripHits);
0128     pixelTrackDeDxHitAss->setValue(i, pixelHits);
0129     stripTrackDeDxHitAss->setValue(i, stripHits);
0130   }
0131 
0132   iEvent.put(std::move(pixelTrackDeDxHitAss), "PixelHits");
0133   iEvent.put(std::move(stripTrackDeDxHitAss), "StripHits");
0134 }
0135 
0136 /*****************************************************************************/
0137 void DeDxHitCalibrator::processHitInfo(const reco::DeDxHitInfo& info,
0138                                        const float& trackMomentum,
0139                                        reco::DeDxHitCollection& pixelHits,
0140                                        reco::DeDxHitCollection& stripHits) {
0141   // Energy loss parametrization from https://doi.org/10.1016/j.nima.2012.06.064
0142   static constexpr double a = 0.07;     //energy loss factor in silicon
0143   static constexpr double l0 = 450e-4;  //reference path length in um
0144 
0145   for (size_t i = 0; i < info.size(); i++) {
0146     // Require hits to be complete and compatible
0147     const auto& type = info.type(i);
0148     if (!(type & (1 << reco::DeDxHitInfo::Complete)) || !(type & (1 << reco::DeDxHitInfo::Compatible)))
0149       continue;
0150 
0151     // Effective path length
0152     const auto& pl = info.pathlength(i);
0153     const auto pathLength = pl * (1. + a * std::log(pl / l0));
0154     const DetId& detId = info.detId(i);
0155 
0156     // Strip
0157     if (const auto& stripCluster = info.stripCluster(i)) {
0158       const auto& thickness =
0159           dynamic_cast<const StripGeomDetUnit*>(tkGeom_->idToDet(detId))->surface().bounds().thickness();
0160       const auto& det = getDetId(detId, thickness);
0161       const auto& chipId = ChipId(detId, stripCluster->barycenter() / sistrip::STRIPS_PER_APV);
0162       // Fill measured strip deposits
0163       const auto& thr = dedxCalib_->thr()[det];
0164       std::vector<std::pair<double, int> > b;
0165       b.reserve(stripCluster->amplitudes().size() + 2);
0166       b.emplace_back(thr, kIsBelow);
0167       for (const auto& adc : stripCluster->amplitudes()) {
0168         if (adc > 253)
0169           b.emplace_back(254., kIsOver);
0170         else
0171           b.emplace_back(adc + 0.5, kIsNormal);
0172       }
0173       b.emplace_back(thr, kIsBelow);
0174       // Fit
0175       const auto& result = fitStripCluster(b, dedxCalib_->alpha()[det], 1. / dedxCalib_->sigma()[det]);
0176       // ADC -> e- -> MeV
0177       const auto& energy = correctEnergy(result.first * sistrip::MeVperADCStrip, chipId);
0178       const auto& esigma = correctEnergy(result.second * sistrip::MeVperADCStrip, chipId);
0179       // Compute charge
0180       const auto& charge = pathLength != 0. ? energy / pathLength : 0.;
0181       stripHits.emplace_back(charge, trackMomentum, pathLength, 0, esigma);
0182     }
0183 
0184     // Pixel
0185     if (const auto& pixelCluster = info.pixelCluster(i)) {
0186       const auto& thickness =
0187           dynamic_cast<const PixelGeomDetUnit*>(tkGeom_->idToDet(detId))->surface().bounds().thickness();
0188       const auto& det = getDetId(detId, thickness);
0189       const auto& chipId =
0190           ChipId(detId, (int(pixelCluster->x() / ROCSizeInX) << 3) + int(pixelCluster->y() / ROCSizeInY));
0191       // Collect adc
0192       double delta(0);
0193       bool isSaturated(false);
0194       for (size_t j = 0; j < pixelCluster->pixelADC().size(); j++) {
0195         const auto& elec = pixelCluster->pixelADC()[j];
0196         delta += elec * MeVPerElectron_;
0197         if (isSaturated)
0198           continue;
0199         // ( pos , (x , y) )
0200         const auto& row = pixelCluster->minPixelRow() + pixelCluster->pixelOffset()[2 * j];
0201         const auto& col = pixelCluster->minPixelCol() + pixelCluster->pixelOffset()[2 * j + 1];
0202         // Go back to adc
0203         const auto& DBgain = pixelCalib_.getGain(detId, col, row);
0204         const auto& DBpedestal = pixelCalib_.getPedestal(detId, col, row);
0205         if (elec == std::numeric_limits<uint16_t>::max())
0206           isSaturated = true;
0207         else if (DBgain > 0.) {
0208           double vcal;
0209           const auto& theLayer = (detId.subdetId() == 1) ? tkTopo_->pxbLayer(detId) : 0;
0210           if (theLayer == 1)
0211             vcal = (elec - VCaltoElectronOffset_L1_) / VCaltoElectronGain_L1_;
0212           else
0213             vcal = (elec - VCaltoElectronOffset_) / VCaltoElectronGain_;
0214           const auto adc = std::round(DBpedestal + vcal / DBgain);
0215           if (adc > pixelSaturationThr_)
0216             isSaturated = true;
0217         }
0218       }
0219       // Compute energy
0220       const auto& energy = correctEnergy(delta, chipId);
0221       if (det != PXF && energy <= 50e-3)
0222         continue;
0223       // Estimate sigma for cluster
0224       const unsigned char& nChannels = pixelCluster->pixelADC().size();
0225       const auto esigma = 10e-3 * std::sqrt(nChannels);
0226       // Compute charge
0227       const auto& charge = pathLength != 0. ? energy / pathLength : 0.;
0228       pixelHits.emplace_back(charge, trackMomentum, pathLength, isSaturated, esigma);
0229     }
0230   }
0231 }
0232 
0233 /*****************************************************************************/
0234 int DeDxHitCalibrator::getDetId(const DetId& id, const float& thickness) {
0235   const auto& subdet = id.subdetId() - 1;
0236   if (subdet == TECThin && thickness > 400e-4)
0237     return TECThick;
0238   return subdet;
0239 }
0240 
0241 /*****************************************************************************/
0242 float DeDxHitCalibrator::correctEnergy(const float& energy, const ChipId& detId) {
0243   const auto& g = dedxCalib_->gain().find(detId);
0244   if (applyGain_ && g != dedxCalib_->gain().end())
0245     return energy * g->second;
0246   return energy;
0247 }
0248 
0249 /*****************************************************************************/
0250 double DeDxHitCalibrator::getChi2(const std::vector<double>& x,
0251                                   const std::vector<std::pair<double, int> >& b,
0252                                   const double& coupling,
0253                                   const double& iSigma) {
0254   const auto& npar = b.size();
0255   std::vector<double> y(npar, 0.);
0256   for (size_t i = 0; i < npar; i++) {
0257     const auto dx = coupling * x[i];
0258     if (i >= 1)
0259       y[i - 1] += dx;
0260     y[i] += x[i] - 2 * dx;
0261     if (i + 1 < npar)
0262       y[i + 1] += dx;
0263   }
0264   double chi2(0.);
0265   for (size_t i = 0; i < npar; i++) {
0266     auto q = (b[i].first - y[i]) * iSigma;
0267     if (b[i].second == kIsNormal)
0268       chi2 += q * q;
0269     else if (b[i].second == kIsBelow) {
0270       q -= 2;
0271       if (q < 0)
0272         chi2 += 0.5 * q * q;
0273     } else if (b[i].second == kIsOver) {
0274       q += 2;
0275       if (q > 0)
0276         chi2 += 0.5 * q * q;
0277     }
0278     // Penalty for negatives
0279     if (x[i] < 0) {
0280       q = x[i] * iSigma;
0281       chi2 += q * q;
0282     }
0283   }
0284   return chi2;
0285 }
0286 
0287 /*****************************************************************************/
0288 void DeDxHitCalibrator::getAlphaBeta(const std::vector<double>& x,
0289                                      const std::vector<std::pair<double, int> >& b,
0290                                      CLHEP::HepMatrix& alpha,
0291                                      CLHEP::HepVector& beta,
0292                                      const std::vector<bool>& isFix,
0293                                      const double& coupling,
0294                                      const double& iSigma) {
0295   const auto& npar = b.size();
0296   const auto a0 = coupling * iSigma;
0297   const auto a1 = (1 - 2 * coupling) * iSigma;
0298   const auto a00 = a0 * a0;
0299   const auto a11 = a1 * a1;
0300   const auto a01 = a0 * a1;
0301   std::vector<double> y(npar, 0.);
0302   for (size_t i = 0; i < npar; i++)
0303     if (!isFix[i]) {
0304       const auto dx = coupling * x[i];
0305       if (i >= 1)
0306         y[i - 1] += dx;
0307       y[i] += x[i] - 2 * dx;
0308       if (i + 1 < npar)
0309         y[i + 1] += dx;
0310     }
0311   for (size_t i = 0; i < npar; i++) {
0312     auto q = (y[i] - b[i].first) * iSigma;
0313     int f(0);
0314     if (b[i].second == kIsNormal)
0315       f = 2;
0316     else if (b[i].second == kIsBelow) {
0317       q += 2;
0318       if (q > 0)
0319         f = 1;
0320     } else if (b[i].second == kIsOver) {
0321       q -= 2;
0322       if (q < 0)
0323         f = 1;
0324     }
0325     if (f > 0) {
0326       if (i >= 1)
0327         if (!isFix[i - 1]) {
0328           alpha[i - 1][i - 1] += f * a00;
0329           if (!isFix[i]) {
0330             alpha[i - 1][i] += f * a01;
0331             alpha[i][i - 1] += f * a01;
0332           }
0333           beta[i - 1] += f * q * a0;
0334         }
0335       if (!isFix[i]) {
0336         alpha[i][i] += f * a11;
0337         beta[i] += f * q * a1;
0338       }
0339       if (i + 1 < npar)
0340         if (!isFix[i + 1]) {
0341           alpha[i + 1][i + 1] += f * a00;
0342           if (!isFix[i]) {
0343             alpha[i + 1][i] += f * a01;
0344             alpha[i][i + 1] += f * a01;
0345           }
0346           beta[i + 1] += f * q * a0;
0347         }
0348     }
0349     // Penalty for negatives
0350     if (!isFix[i])
0351       if (x[i] < 0) {
0352         alpha[i][i] += 2 * iSigma * iSigma;
0353         q = x[i] * iSigma;
0354         beta[i] += 2 * q * iSigma;
0355       }
0356   }
0357   for (size_t i = 0; i < npar; i++)
0358     if (isFix[i]) {
0359       alpha[i][i] = 1.;
0360       beta[i] = 0.;
0361     }
0362 }
0363 
0364 /*****************************************************************************/
0365 std::pair<double, double> DeDxHitCalibrator::fitStripCluster(const std::vector<std::pair<double, int> >& b,
0366                                                              const double& coupling,
0367                                                              const double& iSigma) {
0368   const auto& npar = b.size();
0369   std::vector<bool> isFix(npar, false);
0370   std::vector<double> x;
0371   x.reserve(b.size());
0372   for (const auto& ib : b)
0373     x.emplace_back(ib.first);
0374 
0375   bool ok = false;
0376   CLHEP::HepMatrix hessian(npar, npar);
0377   int iter = 0;
0378   do {
0379     double diff(0);
0380     auto old = getChi2(x, b, coupling, iSigma);
0381     do {
0382       CLHEP::HepMatrix alpha(npar, npar, 0.);
0383       CLHEP::HepVector beta(npar, 0.);
0384       getAlphaBeta(x, b, alpha, beta, isFix, coupling, iSigma);
0385       const auto& delta = CLHEP::solve(alpha, -beta);
0386       double lambda(1.);
0387       std::vector<double> xn(npar);
0388       for (size_t i = 0; i < npar; i++)
0389         xn[i] = x[i] + lambda * delta[i];
0390       const auto& next = getChi2(xn, b, coupling, iSigma);
0391       diff = old - next;
0392       if (diff > 0)
0393         x = xn;
0394       old = next;
0395       hessian = alpha;
0396       iter++;
0397     } while (diff > 1e-6 && iter < 100);
0398     // Check if we have negatives
0399     const auto& mi = std::min_element(x.begin(), x.end()) - x.begin();
0400     if (x[mi] < 0) {
0401       x[mi] = 0.;
0402       isFix[mi] = true;
0403     } else
0404       ok = true;
0405   } while (!ok && iter < 100);
0406   hessian /= 2.;
0407   int flag;
0408   hessian.invert(flag);
0409 
0410   double var2(0.);
0411   for (size_t i = 0; i < npar; i++)
0412     for (size_t j = 0; j < npar; j++)
0413       if (!isFix[i] && !isFix[j])
0414         var2 += hessian[i][j];
0415   var2 = std::abs(var2);
0416 
0417   double sum(0.);
0418   for (size_t i = 0; i < npar; i++)
0419     if (!isFix[i])
0420       sum += x[i];
0421 
0422   return {sum, std::sqrt(var2)};
0423 }
0424 
0425 //define this as a plug-in
0426 #include "FWCore/Framework/interface/MakerMacros.h"
0427 DEFINE_FWK_MODULE(DeDxHitCalibrator);
[ 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!