00001 #include "SiTrivialInduceChargeOnStrips.h"
00002
00003 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00004 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00005 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
00006 #include <Math/ProbFuncMathCore.h>
00007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
00008 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00009 #include "FWCore/Utilities/interface/Exception.h"
00010
00011 #include <algorithm>
00012 #include <iostream>
00013
00014 const int
00015 SiTrivialInduceChargeOnStrips::
00016 Ntypes = 14;
00017
00018 const std::string
00019 SiTrivialInduceChargeOnStrips::
00020 type[Ntypes] = { "IB1", "IB2","OB1","OB2","W1a","W2a","W3a","W1b","W2b","W3b","W4","W5","W6","W7"};
00021
00022 static
00023 std::vector<std::vector<double> >
00024 fillSignalCoupling(const edm::ParameterSet& conf, int nTypes, const std::string* typeArray) {
00025 std::vector<std::vector<double> > signalCoupling;
00026 signalCoupling.reserve(nTypes);
00027 std::string mode = conf.getParameter<bool>("APVpeakmode") ? "Peak" : "Dec";
00028 for(int i=0; i<nTypes; ++i) {
00029 signalCoupling.push_back(conf.getParameter<std::vector<double> >("CouplingConstant"+mode+typeArray[i]));
00030 }
00031 return signalCoupling;
00032 }
00033
00034 inline unsigned int
00035 SiTrivialInduceChargeOnStrips::
00036 indexOf(const std::string& t) { return std::find( type, type + Ntypes, t) - type;}
00037
00038 inline unsigned int
00039 SiTrivialInduceChargeOnStrips::
00040 typeOf(const StripGeomDetUnit& det, const TrackerTopology *tTopo) {
00041 DetId id = det.geographicalId();
00042 switch (det.specificType().subDetector()) {
00043 case GeomDetEnumerators::TIB: {return (tTopo->tibLayer(id) < 3) ? indexOf("IB1") : indexOf("IB2");}
00044 case GeomDetEnumerators::TOB: {return (tTopo->tobLayer(id) > 4) ? indexOf("OB1") : indexOf("OB2");}
00045 case GeomDetEnumerators::TID: {return indexOf("W1a") -1 + tTopo->tidRing(id);}
00046 case GeomDetEnumerators::TEC: {return indexOf("W1b") -1 + tTopo->tecRing(id);}
00047 default: throw cms::Exception("Invalid subdetector") << id();
00048 }
00049 }
00050
00051 SiTrivialInduceChargeOnStrips::
00052 SiTrivialInduceChargeOnStrips(const edm::ParameterSet& conf,double g)
00053 : signalCoupling(fillSignalCoupling(conf, Ntypes, type)), Nsigma(3.), geVperElectron(g) {
00054 }
00055
00056 void
00057 SiTrivialInduceChargeOnStrips::
00058 induce(const SiChargeCollectionDrifter::collection_type& collection_points,
00059 const StripGeomDetUnit& det,
00060 std::vector<double>& localAmplitudes,
00061 size_t& recordMinAffectedStrip,
00062 size_t& recordMaxAffectedStrip,
00063 const TrackerTopology *tTopo) const {
00064
00065 const std::vector<double>& coupling = signalCoupling.at(typeOf(det,tTopo));
00066 const StripTopology& topology = dynamic_cast<const StripTopology&>(det.specificTopology());
00067 size_t Nstrips = topology.nstrips();
00068
00069 for (SiChargeCollectionDrifter::collection_type::const_iterator
00070 signalpoint = collection_points.begin(); signalpoint != collection_points.end(); signalpoint++ ) {
00071
00072
00073 double chargePosition = topology.strip(signalpoint->position());
00074 double chargeSpread = signalpoint->sigma() / topology.localPitch(signalpoint->position());
00075
00076 size_t fromStrip = size_t(std::max( 0, int(std::floor( chargePosition - Nsigma*chargeSpread))));
00077 size_t untilStrip = size_t(std::min( Nstrips, size_t(std::ceil( chargePosition + Nsigma*chargeSpread) )));
00078 for (size_t strip = fromStrip; strip < untilStrip; strip++) {
00079
00080 double chargeDepositedOnStrip = chargeDeposited( strip, Nstrips, signalpoint->amplitude(), chargeSpread, chargePosition);
00081
00082 size_t affectedFromStrip = size_t(std::max( 0, int(strip - coupling.size() + 1)));
00083 size_t affectedUntilStrip = size_t(std::min( Nstrips, strip + coupling.size()) );
00084 for (size_t affectedStrip = affectedFromStrip; affectedStrip < affectedUntilStrip; affectedStrip++) {
00085 localAmplitudes.at( affectedStrip ) += chargeDepositedOnStrip * coupling.at(abs( affectedStrip - strip )) ;
00086 }
00087
00088 if( affectedFromStrip < recordMinAffectedStrip ) recordMinAffectedStrip = affectedFromStrip;
00089 if( affectedUntilStrip > recordMaxAffectedStrip ) recordMaxAffectedStrip = affectedUntilStrip;
00090 }
00091 }
00092 return;
00093 }
00094
00095 inline double
00096 SiTrivialInduceChargeOnStrips::
00097 chargeDeposited(size_t strip, size_t Nstrips, double amplitude, double chargeSpread, double chargePosition) const {
00098 double integralUpToStrip = (strip == 0) ? 0. : ( ROOT::Math::normal_cdf( strip, chargeSpread, chargePosition) );
00099 double integralUpToNext = (strip+1 == Nstrips) ? 1. : ( ROOT::Math::normal_cdf( strip+1, chargeSpread, chargePosition) );
00100 double percentOfSignal = integralUpToNext - integralUpToStrip;
00101
00102 return percentOfSignal * amplitude / geVperElectron;
00103 }