SiStripDigitizer/plugins/SiTrivialInduceChargeOnStrips.cc

0001 #include "SiTrivialInduceChargeOnStrips.h"
0002
0003 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0004 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0005 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0006 #include "DataFormats/Math/interface/approx_erf.h"
0007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0008 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010
0011 #include <algorithm>
0012 #include <iostream>
0013
0014 namespace {
0015   struct Count {
0016     Count() {}
0017 #ifdef SISTRIP_COUNT
0018     // note: this code is not thread safe, counts will be inaccurate if run with multiple threads
0019     double ncall = 0;
0020     double ndep = 0, ndep2 = 0, maxdep = 0;
0021     double nstr = 0, nstr2 = 0;
0022     double ncv = 0, nval = 0, nval2 = 0, maxv = 0;
0023     double dzero = 0;
0024     void dep(double d) {
0025       ncall++;
0026       ndep += d;
0027       ndep2 += d * d;
0028       maxdep = std::max(d, maxdep);
0029     }
0030     void str(double d) {
0031       nstr += d;
0032       nstr2 += d * d;
0033     }
0034     void val(double d) {
0035       ncv++;
0036       nval += d;
0037       nval2 += d * d;
0038       maxv = std::max(d, maxv);
0039     }
0040     void zero() { dzero++; }
0041     ~Count() {
0042       std::cout << "deposits " << ncall << " " << maxdep << " " << ndep / ncall << " "
0043                 << std::sqrt(ndep2 * ncall - ndep * ndep) / ncall << std::endl;
0044       std::cout << "zeros " << dzero << std::endl;
0045       std::cout << "strips  " << nstr / ndep << " " << std::sqrt(nstr2 * ndep - nstr * nstr) / ndep << std::endl;
0046       std::cout << "vaules  " << ncv << " " << maxv << " " << nval / ncv << " "
0047                 << std::sqrt(nval2 * ncv - nval * nval) / ncv << std::endl;
0048     }
0049   };
0050   Count count;
0051 #else
0052     void dep(double) const {}
0053     void str(double) const {}
0054     void val(double) const {}
0055     void zero() const {}
0056   };
0057   const Count count;
0058 #endif
0059
0060 }  // namespace
0061
0062 namespace {
0063   constexpr int Ntypes = 14;
0064   //                                   0     1      2     3     4    5      6     7
0065   const std::string type[Ntypes] = {
0066       "IB1", "IB2", "OB1", "OB2", "W1a", "W2a", "W3a", "W1b", "W2b", "W3b", "W4", "W5", "W6", "W7"};
0067   enum { indexOfIB1 = 0, indexOfIB2 = 1, indexOfOB1 = 2, indexOfOB2 = 3, indexOfW1a = 4, indexOfW1b = 7 };
0068
0069   inline std::vector<std::vector<float> > fillSignalCoupling(const edm::ParameterSet& conf,
0070                                                              int nTypes,
0071                                                              const std::string* typeArray) {
0072     std::vector<std::vector<float> > signalCoupling;
0073     signalCoupling.reserve(nTypes);
0074     std::string mode = conf.getParameter<bool>("APVpeakmode") ? "Peak" : "Dec";
0075     for (int i = 0; i < nTypes; ++i) {
0076       std::string version = "";
0077       if (typeArray[i].find('W') != std::string::npos && mode == "Dec")
0078         version = conf.getParameter<bool>("CouplingConstantsRunIIDecW") ? "RunII" : "";
0079       else if (typeArray[i].find('B') != std::string::npos && mode == "Dec")
0080         version = conf.getParameter<bool>("CouplingConstantsRunIIDecB") ? "RunII" : "";
0081
0082       auto dc = conf.getParameter<std::vector<double> >("CouplingConstant" + version + mode + typeArray[i]);
0083       signalCoupling.emplace_back(dc.begin(), dc.end());
0084     }
0085     return signalCoupling;
0086   }
0087
0088   inline unsigned int typeOf(const StripGeomDetUnit& det, const TrackerTopology* tTopo) {
0089     DetId id = det.geographicalId();
0090     switch (det.specificType().subDetector()) {
0091       case GeomDetEnumerators::TIB: {
0092         return (tTopo->tibLayer(id) < 3) ? indexOfIB1 : indexOfIB2;
0093       }
0094       case GeomDetEnumerators::TOB: {
0095         return (tTopo->tobLayer(id) > 4) ? indexOfOB1 : indexOfOB2;
0096       }
0097       case GeomDetEnumerators::TID: {
0098         return indexOfW1a - 1 + tTopo->tidRing(id);
0099       }  //fragile: relies on ordering of 'type'
0100       case GeomDetEnumerators::TEC: {
0101         return indexOfW1b - 1 + tTopo->tecRing(id);
0102       }  //fragile: relies on ordering of 'type'
0103       default:
0104         throw cms::Exception("Invalid subdetector") << id();
0105     }
0106   }
0107
0108   inline double chargeDeposited(
0109       size_t strip, size_t Nstrips, double amplitude, double chargeSpread, double chargePosition) {
0110     double integralUpToStrip =
0111         (strip == 0) ? 0. : (approx_erf((strip - chargePosition) / chargeSpread / 1.41421356237309515));
0112     double integralUpToNext =
0113         (strip + 1 == Nstrips) ? 1. : (approx_erf((strip + 1 - chargePosition) / chargeSpread / 1.41421356237309515));
0114
0115     return 0.5 * (integralUpToNext - integralUpToStrip) * amplitude;
0116   }
0117
0118 }  // namespace
0119
0120 SiTrivialInduceChargeOnStrips::SiTrivialInduceChargeOnStrips(const edm::ParameterSet& conf, double g)
0121     : signalCoupling(fillSignalCoupling(conf, Ntypes, type)), Nsigma(3.), geVperElectron(g) {}
0122
0123 void SiTrivialInduceChargeOnStrips::induce(const SiChargeCollectionDrifter::collection_type& collection_points,
0124                                            const StripGeomDetUnit& det,
0125                                            std::vector<float>& localAmplitudes,
0126                                            size_t& recordMinAffectedStrip,
0127                                            size_t& recordMaxAffectedStrip,
0128                                            const TrackerTopology* tTopo) const {
0129   induceVector(collection_points, det, localAmplitudes, recordMinAffectedStrip, recordMaxAffectedStrip, tTopo);
0130
0131   /*
0132   auto ominA=recordMinAffectedStrip, omaxA=recordMaxAffectedStrip;
0133   std::vector<float> oampl(localAmplitudes);
0134   induceOriginal(collection_points, det, oampl, ominA, omaxA, tTopo);
0135
0136   //  std::cout << "orig " << ominA << " " << omaxA << " ";
0137   //for (auto a : oampl) std::cout << a << ",";
0138   //std::cout << std::endl;
0139
0140   auto minA=recordMinAffectedStrip, maxA=recordMaxAffectedStrip;
0141   std::vector<float> ampl(localAmplitudes);
0142   induceVector(collection_points, det, ampl, minA, maxA, tTopo);
0143
0144   // std::cout << "vect " << minA << " " << maxA << " ";
0145   //for (auto a :   ampl) std::cout << a << ",";
0146   //std::cout << std::endl;
0147
0148   float diff=0;
0149   for (size_t i=0; i!=ampl.size(); ++i) { diff = std::max(diff,ampl[i]>0 ? std::abs(ampl[i]-oampl[i])/ampl[i] : 0);}
0150   if (diff> 1.e-4) {
0151     std::cout << diff << std::endl;
0152     std::cout << "orig " << ominA << " " << omaxA << " ";
0153 //    for (auto a : oampl) std::cout << a << ",";
0154     std::cout << std::endl;
0155     std::cout << "vect " << minA << " " << maxA << " ";
0156 //    for (auto a : ampl) std::cout << a << ",";
0157     std::cout << std::endl;
0158   }
0159
0160   localAmplitudes.swap(ampl);
0161   recordMinAffectedStrip=minA;
0162   recordMaxAffectedStrip=maxA;
0163   */
0164 }
0165
0166 void SiTrivialInduceChargeOnStrips::induceVector(const SiChargeCollectionDrifter::collection_type& collection_points,
0167                                                  const StripGeomDetUnit& det,
0168                                                  std::vector<float>& localAmplitudes,
0169                                                  size_t& recordMinAffectedStrip,
0170                                                  size_t& recordMaxAffectedStrip,
0171                                                  const TrackerTopology* tTopo) const {
0172   auto const& coupling = signalCoupling[typeOf(det, tTopo)];
0173   const StripTopology& topology = dynamic_cast<const StripTopology&>(det.specificTopology());
0174   const int Nstrips = topology.nstrips();
0175
0176   if (Nstrips == 0)
0177     return;
0178
0179   const int NP = collection_points.size();
0180   if (0 == NP)
0181     return;
0182
0183   constexpr int MaxN = 512;
0184   // if NP too large split...
0185
0186   for (int ip = 0; ip < NP; ip += MaxN) {
0187     auto N = std::min(NP - ip, MaxN);
0188
0189     count.dep(N);
0190     float amplitude[N];
0191     float chargePosition[N];
0192     float chargeSpread[N];
0193     int fromStrip[N];
0194     int nStrip[N];
0195
0196     // load not vectorize
0197     //In strip coordinates:
0198     for (int i = 0; i != N; ++i) {
0199       auto j = ip + i;
0200       if (0 == collection_points[j].amplitude())
0201         count.zero();
0202       chargePosition[i] = topology.strip(collection_points[j].position());
0203       chargeSpread[i] = collection_points[j].sigma() / topology.localPitch(collection_points[j].position());
0204       amplitude[i] = 0.5f * collection_points[j].amplitude() / geVperElectron;
0205     }
0206
0207     // this vectorize
0208     for (int i = 0; i != N; ++i) {
0209       fromStrip[i] = std::max(0, int(std::floor(chargePosition[i] - Nsigma * chargeSpread[i])));
0210       nStrip[i] = std::min(Nstrips, int(std::ceil(chargePosition[i] + Nsigma * chargeSpread[i]))) - fromStrip[i];
0211     }
0212     int tot = 0;
0213     for (int i = 0; i != N; ++i)
0214       tot += nStrip[i];
0215     tot += N;  // add last strip
0216     count.val(tot);
0217     float value[tot];
0218
0219     // assign relative position (lower bound of strip) in value;
0220     int kk = 0;
0221     for (int i = 0; i != N; ++i) {
0222       auto delta = 1.f / (std::sqrt(2.f) * chargeSpread[i]);
0223       auto pos = delta * (float(fromStrip[i]) - chargePosition[i]);
0224
0225       // VI: before value[0] was not defined and value[tot] was filled
0226       //     to fix this the loop below was changed
0227       for (int j = 0; j <= nStrip[i]; ++j) {  /// include last strip
0228         value[kk] = pos + float(j) * delta;
0229         ++kk;
0230       }
0231     }
0232     assert(kk == tot);
0233
0234     // main loop fully vectorized
0235     for (int k = 0; k != tot; ++k)
0236       value[k] = approx_erf(value[k]);
0237
0238     // saturate 0 & NStrips strip to 0 and 1???
0239     kk = 0;
0240     for (int i = 0; i != N; ++i) {
0241       if (0 == fromStrip[i])
0242         value[kk] = 0;
0243       kk += nStrip[i];
0244       if (Nstrips == fromStrip[i] + nStrip[i])
0245         value[kk] = 1.f;
0246       ++kk;
0247     }
0248     assert(kk == tot);
0249
0250     // compute integral over strip (lower bound becomes the value)
0251     for (int k = 0; k != tot - 1; ++k)
0252       value[k] -= value[k + 1];  // this is negative!
0253
0254     float charge[Nstrips];
0255     for (int i = 0; i != Nstrips; ++i)
0256       charge[i] = 0;
0257     kk = 0;
0258     for (int i = 0; i != N; ++i) {
0259       for (int j = 0; j != nStrip[i]; ++j)
0260         charge[fromStrip[i] + j] -= amplitude[i] * value[kk++];
0261       ++kk;  // skip last "strip"
0262     }
0263     assert(kk == tot);
0264
0265     /// do crosstalk... (can be done better, most probably not worth)
0266     int minA = recordMinAffectedStrip, maxA = recordMaxAffectedStrip;
0267     int sc = coupling.size();
0268     for (int i = 0; i != Nstrips; ++i) {
0269       int strip = i;
0270       if (0 == charge[i])
0271         continue;
0272       auto affectedFromStrip = std::max(0, strip - sc + 1);
0273       auto affectedUntilStrip = std::min(Nstrips, strip + sc);
0274       for (auto affectedStrip = affectedFromStrip; affectedStrip < affectedUntilStrip; ++affectedStrip)
0275         localAmplitudes[affectedStrip] += charge[i] * coupling[std::abs(affectedStrip - strip)];
0276
0277       if (affectedFromStrip < minA)
0278         minA = affectedFromStrip;
0279       if (affectedUntilStrip > maxA)
0280         maxA = affectedUntilStrip;
0281     }
0282     recordMinAffectedStrip = minA;
0283     recordMaxAffectedStrip = maxA;
0284   }  // end loop ip
0285 }
0286
0287 void SiTrivialInduceChargeOnStrips::induceOriginal(const SiChargeCollectionDrifter::collection_type& collection_points,
0288                                                    const StripGeomDetUnit& det,
0289                                                    std::vector<float>& localAmplitudes,
0290                                                    size_t& recordMinAffectedStrip,
0291                                                    size_t& recordMaxAffectedStrip,
0292                                                    const TrackerTopology* tTopo) const {
0293   auto const& coupling = signalCoupling[typeOf(det, tTopo)];
0294   const StripTopology& topology = dynamic_cast<const StripTopology&>(det.specificTopology());
0295   size_t Nstrips = topology.nstrips();
0296
0297   if (!collection_points.empty())
0298     count.dep(collection_points.size());
0299
0300   for (auto signalpoint = collection_points.begin(); signalpoint != collection_points.end(); signalpoint++) {
0301     //In strip coordinates:
0302     double chargePosition = topology.strip(signalpoint->position());
0303     double chargeSpread = signalpoint->sigma() / topology.localPitch(signalpoint->position());
0304
0305     size_t fromStrip = size_t(std::max(0, int(std::floor(chargePosition - Nsigma * chargeSpread))));
0306     size_t untilStrip = size_t(std::min(Nstrips, size_t(std::ceil(chargePosition + Nsigma * chargeSpread))));
0307
0308     count.str(std::max(0, int(untilStrip) - int(fromStrip)));
0309     for (size_t strip = fromStrip; strip < untilStrip; strip++) {
0310       double chargeDepositedOnStrip =
0311           chargeDeposited(strip, Nstrips, signalpoint->amplitude() / geVperElectron, chargeSpread, chargePosition);
0312
0313       size_t affectedFromStrip = size_t(std::max(0, int(strip - coupling.size() + 1)));
0314       size_t affectedUntilStrip = size_t(std::min(Nstrips, strip + coupling.size()));
0315       for (size_t affectedStrip = affectedFromStrip; affectedStrip < affectedUntilStrip; affectedStrip++) {
0316         localAmplitudes.at(affectedStrip) +=
0317             chargeDepositedOnStrip * coupling.at((size_t)std::abs((int)affectedStrip - (int)strip));
0318       }
0319
0320       if (affectedFromStrip < recordMinAffectedStrip)
0321         recordMinAffectedStrip = affectedFromStrip;
0322       if (affectedUntilStrip > recordMaxAffectedStrip)
0323         recordMaxAffectedStrip = affectedUntilStrip;
0324     }
0325   }
0326 }