#include <SiTrivialInduceChargeOnStrips.h>

Inheritance diagram for SiTrivialInduceChargeOnStrips:

Public Member Functions
void	induce (const SiChargeCollectionDrifter::collection_type &collection_points, const StripGeomDetUnit &det, std::vector< float > &localAmplitudes, size_t &recordMinAffectedStrip, size_t &recordMaxAffectedStrip, const TrackerTopology *tTopo) const

	SiTrivialInduceChargeOnStrips (const edm::ParameterSet &conf, double g)

virtual	~SiTrivialInduceChargeOnStrips ()

Public Member Functions inherited from SiInduceChargeOnStrips
virtual	~SiInduceChargeOnStrips ()

Private Member Functions
void	induceOriginal (const SiChargeCollectionDrifter::collection_type &collection_points, const StripGeomDetUnit &det, std::vector< float > &localAmplitudes, size_t &recordMinAffectedStrip, size_t &recordMaxAffectedStrip, const TrackerTopology *tTopo) const

void	induceVector (const SiChargeCollectionDrifter::collection_type &collection_points, const StripGeomDetUnit &det, std::vector< float > &localAmplitudes, size_t &recordMinAffectedStrip, size_t &recordMaxAffectedStrip, const TrackerTopology *tTopo) const

Private Attributes
const float	geVperElectron

const float	Nsigma

const std::vector< std::vector < float > >	signalCoupling

Detailed Description

Definition at line 9 of file SiTrivialInduceChargeOnStrips.h.

Constructor & Destructor Documentation

SiTrivialInduceChargeOnStrips::SiTrivialInduceChargeOnStrips	(	const edm::ParameterSet &	conf,
		double	g
	)

Definition at line 90 of file SiTrivialInduceChargeOnStrips.cc.

91 : signalCoupling(fillSignalCoupling(conf, Ntypes, type)), Nsigma(3.), geVperElectron(g) {

92 }

type

Definition: HCALResponse.h:21

SiTrivialInduceChargeOnStrips::signalCoupling

const std::vector< std::vector< float > > signalCoupling

Definition: SiTrivialInduceChargeOnStrips.h:38

g

The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g

Definition: Activities.doc:4

SiTrivialInduceChargeOnStrips::Nsigma

const float Nsigma

Definition: SiTrivialInduceChargeOnStrips.h:40

SiTrivialInduceChargeOnStrips::geVperElectron

const float geVperElectron

Definition: SiTrivialInduceChargeOnStrips.h:41

virtual SiTrivialInduceChargeOnStrips::~SiTrivialInduceChargeOnStrips ( )

inlinevirtual

Definition at line 12 of file SiTrivialInduceChargeOnStrips.h.

12 {}

Member Function Documentation

void SiTrivialInduceChargeOnStrips::induce	(	const SiChargeCollectionDrifter::collection_type &	collection_points,
		const StripGeomDetUnit &	det,
		std::vector< float > &	localAmplitudes,
		size_t &	recordMinAffectedStrip,
		size_t &	recordMaxAffectedStrip,
		const TrackerTopology *	tTopo
	)		const

virtual

Implements SiInduceChargeOnStrips.

Definition at line 96 of file SiTrivialInduceChargeOnStrips.cc.

References induceVector().

                                            {
 
 
    induceVector(collection_points, det, localAmplitudes, recordMinAffectedStrip, recordMaxAffectedStrip, tTopo);
 
   /*
   auto ominA=recordMinAffectedStrip, omaxA=recordMaxAffectedStrip;
   std::vector<float> oampl(localAmplitudes);    
   induceOriginal(collection_points, det, oampl, ominA, omaxA, tTopo);
 
   //  std::cout << "orig " << ominA << " " << omaxA << " ";
   //for (auto a : oampl) std::cout << a << ",";
   //std::cout << std::endl;
 
   auto minA=recordMinAffectedStrip, maxA=recordMaxAffectedStrip;
   std::vector<float> ampl(localAmplitudes);
   induceVector(collection_points, det, ampl, minA, maxA, tTopo);
 
   // std::cout << "vect " << minA << " " << maxA << " ";          
   //for (auto a :   ampl) std::cout << a << ",";
   //std::cout << std::endl;
  
   float diff=0;
   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);}
   if (diff> 1.e-4) {
     std::cout << diff << std::endl;
     std::cout << "orig " << ominA << " " << omaxA << " ";
 //    for (auto a : oampl) std::cout << a << ",";
     std::cout << std::endl;
     std::cout << "vect " << minA << " " << maxA << " ";
 //    for (auto a : ampl) std::cout << a << ",";
     std::cout << std::endl;
   }
 
   localAmplitudes.swap(ampl);
   recordMinAffectedStrip=minA;
   recordMaxAffectedStrip=maxA;
   */
 
 }

void SiTrivialInduceChargeOnStrips::induceOriginal	(	const SiChargeCollectionDrifter::collection_type &	collection_points,
		const StripGeomDetUnit &	det,
		std::vector< float > &	localAmplitudes,
		size_t &	recordMinAffectedStrip,
		size_t &	recordMaxAffectedStrip,
		const TrackerTopology *	tTopo
	)		const

private

Definition at line 261 of file SiTrivialInduceChargeOnStrips.cc.

References funct::abs(), KineDebug3::count(), geVperElectron, StripTopology::localPitch(), bookConverter::max, min(), Nsigma, StripTopology::nstrips(), StripGeomDetUnit::specificTopology(), StripTopology::strip(), and ecaldqm::topology().

                                                {
 
 
   auto const & coupling = signalCoupling[typeOf(det,tTopo)];
   const StripTopology& topology = dynamic_cast<const StripTopology&>(det.specificTopology());
   size_t Nstrips =  topology.nstrips();
 
   if (!collection_points.empty()) count.dep(collection_points.size());
 
   for (auto signalpoint = collection_points.begin();  signalpoint != collection_points.end();  signalpoint++ ) {
     
     //In strip coordinates:
     double chargePosition = topology.strip(signalpoint->position());
     double chargeSpread = signalpoint->sigma() / topology.localPitch(signalpoint->position());
     
     size_t fromStrip  = size_t(std::max( 0,          int(std::floor( chargePosition - Nsigma*chargeSpread))));
     size_t untilStrip = size_t(std::min( Nstrips, size_t(std::ceil( chargePosition + Nsigma*chargeSpread) )));
 
     count.str(std::max(0,int(untilStrip)-int(fromStrip)));
     for (size_t strip = fromStrip;  strip < untilStrip; strip++) {
 
       double chargeDepositedOnStrip = chargeDeposited( strip, Nstrips, signalpoint->amplitude() / geVperElectron, chargeSpread, chargePosition);
 
       size_t affectedFromStrip  = size_t(std::max( 0, int(strip - coupling.size() + 1)));
       size_t affectedUntilStrip = size_t(std::min( Nstrips, strip + coupling.size())   );  
       for (size_t affectedStrip = affectedFromStrip;  affectedStrip < affectedUntilStrip;  affectedStrip++) {
     localAmplitudes.at( affectedStrip ) += chargeDepositedOnStrip * coupling.at(abs( affectedStrip - strip )) ;
       }
 
       if( affectedFromStrip  < recordMinAffectedStrip ) recordMinAffectedStrip = affectedFromStrip;
       if( affectedUntilStrip > recordMaxAffectedStrip ) recordMaxAffectedStrip = affectedUntilStrip;
     }
   }
 
 }

void SiTrivialInduceChargeOnStrips::induceVector	(	const SiChargeCollectionDrifter::collection_type &	collection_points,
		const StripGeomDetUnit &	det,
		std::vector< float > &	localAmplitudes,
		size_t &	recordMinAffectedStrip,
		size_t &	recordMaxAffectedStrip,
		const TrackerTopology *	tTopo
	)		const

private

include last strip

do crosstalk... (can be done better, most probably not worth)

Definition at line 144 of file SiTrivialInduceChargeOnStrips.cc.

References funct::abs(), approx_erf(), assert(), RecoTauCleanerPlugins::charge, constexpr, KineDebug3::count(), delta, f, geVperElectron, i, j, relval_2017::k, GetRecoTauVFromDQM_MC_cff::kk, StripTopology::localPitch(), bookConverter::max, HLT_25ns10e33_v2_cff::MaxN, min(), N, Nsigma, StripTopology::nstrips(), position, StripGeomDetUnit::specificTopology(), mathSSE::sqrt(), StripTopology::strip(), ecaldqm::topology(), and relativeConstraints::value.

Referenced by induce().

                                                {
 
 
   auto const & coupling = signalCoupling[typeOf(det,tTopo)];
   const StripTopology& topology = dynamic_cast<const StripTopology&>(det.specificTopology());
   const int Nstrips =  topology.nstrips();
 
   if (Nstrips == 0) return;
 
   const int NP = collection_points.size();
   if(0==NP) return;
 
   constexpr int MaxN = 512;
   // if NP too large split...
 
   for (int ip=0; ip<NP; ip+=MaxN) {
     auto N = std::min(NP-ip,MaxN);
 
     count.dep(N);
     float amplitude[N];
     float chargePosition[N];
     float chargeSpread[N];
     int fromStrip[N];
     int nStrip[N];
 
     // load not vectorize
     //In strip coordinates:
     for (int i=0; i!=N;++i) {
       auto j = ip+i;
       if (0==collection_points[j].amplitude()) count.zero();
       chargePosition[i]=topology.strip(collection_points[j].position());
       chargeSpread[i]= collection_points[j].sigma() / topology.localPitch(collection_points[j].position());
       amplitude[i]=0.5f*collection_points[j].amplitude() / geVperElectron;
     }
     
     // this vectorize
     for (int i=0; i!=N;++i) {
       fromStrip[i]  = std::max( 0,  int(std::floor( chargePosition[i] - Nsigma*chargeSpread[i])) );
       nStrip[i] = std::min( Nstrips, int(std::ceil( chargePosition[i] + Nsigma*chargeSpread[i])) ) - fromStrip[i];
     }
     int tot=0;
     for (int i=0; i!=N;++i) tot += nStrip[i];
     tot+=N; // add last strip 
     count.val(tot);
     float value[tot];
     
     // assign relative position (lower bound of strip) in value;
     int kk=0;
     for (int i=0; i!=N;++i) {
       auto delta = 1.f/(std::sqrt(2.f)*chargeSpread[i]);
       auto pos = delta*(float(fromStrip[i])-chargePosition[i]);
 
       // VI: before value[0] was not defined and value[tot] was filled
       //     to fix this the loop below was changed
       for (int j=0;j<=nStrip[i]; ++j) { 
     value[kk] = pos+float(j)*delta;
         ++kk;  
       }
     }
     assert(kk==tot);
     
     // main loop fully vectorized
     for (int k=0;k!=tot; ++k)
       value[k] = approx_erf(value[k]);
     
     // saturate 0 & NStrips strip to 0 and 1???
     kk=0;
     for (int i=0; i!=N;++i) {
       if (0 == fromStrip[i])  value[kk]=0;
       kk+=nStrip[i];
       if (Nstrips == fromStrip[i]+nStrip[i]) value[kk]=1.f;
       ++kk;
     }
     assert(kk==tot);
     
     // compute integral over strip (lower bound becomes the value)
     for (int k=0;k!=tot-1; ++k)
       value[k]-=value[k+1];  // this is negative!
     
     
     float charge[Nstrips]; for (int i=0;i!=Nstrips; ++i) charge[i]=0;
     kk=0;
     for (int i=0; i!=N;++i){ 
       for (int j=0;j!=nStrip[i]; ++j)
     charge[fromStrip[i]+j]-= amplitude[i]*value[kk++];
       ++kk; // skip last "strip"
     }
     assert(kk==tot);
     
     
     int minA=recordMinAffectedStrip, maxA=recordMaxAffectedStrip;
     int sc = coupling.size();
     for (int i=0;i!=Nstrips; ++i) {
       int strip = i;
       if (0==charge[i]) continue;
       auto affectedFromStrip  = std::max( 0, strip - sc + 1);
       auto affectedUntilStrip = std::min(Nstrips, strip + sc);  
       for (auto affectedStrip=affectedFromStrip;  affectedStrip < affectedUntilStrip;  ++affectedStrip)
     localAmplitudes[affectedStrip] += charge[i] * coupling[std::abs(affectedStrip - strip)] ;
       
       if( affectedFromStrip  < minA ) minA = affectedFromStrip;
       if( affectedUntilStrip > maxA ) maxA = affectedUntilStrip;
     }
     recordMinAffectedStrip=minA;
     recordMaxAffectedStrip=maxA;
   }  // end loop ip
 
 }