9 #include <gsl/gsl_sf_erf.h> 23 signalCoupling.clear();
29 int stripLeft, stripRight, stripLeftW, stripRightW;
30 double upperBound, lowerBound, upperBoundW, lowerBoundW;
35 std::map<int, float, less<int>> x, y;
37 for (CDrifterFP420::collection_type::const_iterator sp = _collection_points.begin(); sp != _collection_points.end();
39 float chargePositionW = -1.;
40 float chargePosition = -1.;
43 G4ThreeVector Position3D = (*sp).position();
46 std::cout <<
" =============================*InduceChargeFP420:induce:" 48 << Position3D << std::endl;
49 std::cout <<
" localPitch= " << localPitch << std::endl;
50 std::cout <<
" xytype= " << xytype <<
" numStrips= " << numStrips << std::endl;
59 std::cout <<
" *InduceChargeFP420:Z slice outside the plate: Position3D= " << Position3D << std::endl;
63 std::cout <<
" *InduceChargeFP420:XY slice outside the plate: Position3D= " << Position3D << std::endl;
70 chargePosition = 0.5 * numStrips + Position3D.x() / localPitch;
72 chargePositionW = 0.5 * numStripsW + Position3D.y() / localPitchW;
76 else if (xytype == 2) {
78 chargePosition = 0.5 * numStrips + Position3D.y() / localPitch;
80 chargePositionW = 0.5 * numStripsW + Position3D.x() / localPitchW;
83 std::cout <<
"**** InduceChargeFP420: !!! ERROR: you have not to be " 85 << xytype << std::endl;
91 std::cout <<
"===================================**** " 92 "InduceChargeFP420: xytype= " 93 << xytype << std::endl;
94 std::cout <<
" chargePositionW= " << chargePositionW <<
" chargePosition= " << chargePosition << std::endl;
95 std::cout <<
"Position3D= " << Position3D << std::endl;
99 float chargeSpread = (*sp).sigma() / localPitch;
100 float chargeSpreadW = (*sp).sigma() / localPitchW;
104 stripLeft =
int(chargePosition - clusterWidth * chargeSpread);
105 stripRight =
int(chargePosition + clusterWidth * chargeSpread);
106 stripLeft = (0 < stripLeft ? stripLeft : 0);
107 stripRight = (numStrips > stripRight ? stripRight : numStrips - 1);
109 stripLeftW =
int(chargePositionW - clusterWidth * chargeSpreadW);
110 stripRightW =
int(chargePositionW + clusterWidth * chargeSpreadW);
111 stripLeftW = (0 < stripLeftW ? stripLeftW : 0);
112 stripRightW = (numStripsW > stripRightW ? stripRightW : numStripsW - 1);
115 std::cout <<
" Position3D = " << Position3D <<
"amplitude=" << (*sp).amplitude() << std::endl;
116 std::cout <<
" MaxChargeSpread= " << clusterWidth * chargeSpread <<
" sigma= " << (*sp).sigma() << std::endl;
117 std::cout <<
"*** numStrips= " << numStrips <<
" localPitch= " << localPitch <<
"xytype=" << xytype << std::endl;
118 std::cout <<
" chargePosition= " << chargePosition <<
" chargeSpread= " << chargeSpread << std::endl;
119 std::cout <<
" stripLeft= " << stripLeft <<
" stripRight= " << stripRight << std::endl;
120 std::cout <<
" chargePositionW= " << chargePositionW <<
" chargeSpreadW= " << chargeSpreadW << std::endl;
121 std::cout <<
" stripLeftW= " << stripLeftW <<
" stripRightW= " << stripRightW << std::endl;
125 for (
int i = stripLeft;
i <= stripRight;
i++) {
132 int status = gsl_sf_erf_Q_e((
i - chargePosition) / chargeSpread, &result);
134 std::cerr <<
"InduceChargeFP420::could not compute gaussian tail " 135 "probability for the threshold chosen" 137 lowerBound = 1. - result.val;
139 std::cout <<
"go to left: i= " <<
i <<
"lowerBound=" << lowerBound << std::endl;
144 if (
i == numStrips - 1)
148 int status = gsl_sf_erf_Q_e((
i - chargePosition + 1) / chargeSpread, &result);
150 std::cerr <<
"InduceChargeFP420::could not compute gaussian tail " 151 "probability for the threshold chosen" 153 upperBound = 1. - result.val;
156 std::cout <<
"go to right: i= " <<
i <<
"upperBound=" << upperBound << std::endl;
160 double totalIntegrationRange = upperBound - lowerBound;
161 x[
i] = totalIntegrationRange;
163 if (totalIntegrationRange <= 0.)
164 std::cout <<
" upperBound= " << upperBound <<
" lowerBound= " << lowerBound << std::endl;
165 if (verbosity == -30) {
166 std::cout <<
" *InduceChargeFP420:====================================X i = " <<
i << std::endl;
167 std::cout <<
" upperBound= " << upperBound <<
" lowerBound= " << lowerBound << std::endl;
168 std::cout <<
" totalIntegrationRange= " << totalIntegrationRange << std::endl;
169 std::cout <<
" *InduceChargeFP420:==================================== " << std::endl;
175 for (
int i = stripLeftW;
i <= stripRightW;
i++) {
182 int status = gsl_sf_erf_Q_e((
i - chargePositionW) / chargeSpreadW, &result);
184 std::cerr <<
"InduceChargeFP420::W could not compute gaussian tail " 185 "probability for the threshold chosen" 187 lowerBoundW = 1. - result.val;
190 std::cout <<
"go to left: i= " <<
i <<
"lowerBoundW=" << lowerBoundW << std::endl;
195 if (
i == numStripsW - 1)
199 int status = gsl_sf_erf_Q_e((
i - chargePositionW + 1) / chargeSpreadW, &result);
201 std::cerr <<
"InduceChargeFP420::W could not compute gaussian tail " 202 "probability for the threshold chosen" 204 upperBoundW = 1. - result.val;
207 std::cout <<
"go to right: i= " <<
i <<
"upperBoundW=" << upperBoundW << std::endl;
211 double totalIntegrationRange = upperBoundW - lowerBoundW;
212 y[
i] = totalIntegrationRange;
214 if (totalIntegrationRange <= 0.)
215 std::cout <<
" upperBoundW= " << upperBoundW <<
" lowerBoundW= " << lowerBoundW << std::endl;
217 if (verbosity == -30) {
218 std::cout <<
" *InduceChargeFP420:====================================XW i= " <<
i << std::endl;
219 std::cout <<
" upperBoundW= " << upperBoundW <<
" lowerBoundW= " << lowerBoundW << std::endl;
220 std::cout <<
" totalIntegrationRange= " << totalIntegrationRange << std::endl;
221 std::cout <<
" *InduceChargeFP420:==================================== " << std::endl;
226 int nSignalCoupling = signalCoupling.size();
240 "*************************************************************** " 242 std::cout <<
" numStripsW= " << numStripsW <<
" numStrips= " << numStrips << std::endl;
243 std::cout <<
" nSignalCoupling= " << nSignalCoupling <<
" xytype= " << xytype << std::endl;
244 std::cout <<
" stripLeftW= " << stripLeftW <<
" stripRightW= " << stripRightW << std::endl;
245 std::cout <<
" stripLeft= " << stripLeft <<
" stripRight= " << stripRight << std::endl;
248 for (
int iy = stripLeftW; iy <= stripRightW; iy++) {
249 for (
int ix = stripLeft; ix <= stripRight; ix++) {
250 for (
int k = -nSignalCoupling + 1;
k <= nSignalCoupling - 1;
k++) {
251 if (ix +
k >= 0 && ix +
k < numStrips) {
252 float ChargeFraction = signalCoupling[
abs(
k)] * (x[ix] * y[iy] / geVperElectron) * (*sp).amplitude();
253 if (ChargeFraction > 0.) {
255 int chan = iy * numStrips + (ix +
k);
263 if (
k == 0 && xytype == 2) {
266 << chan << std::endl;
267 std::cout <<
"ix= " << ix <<
"iy= " << iy <<
"k= " <<
k <<
"ChargeFraction= " << ChargeFraction
269 std::cout <<
"hit_signal[chan]= " << hit_signal[
chan] <<
"geVperElectron= " << geVperElectron
271 std::cout <<
"signalCoupling[abs(k)]= " << signalCoupling[
abs(
k)] <<
"x[ix]= " << x[ix]
272 <<
"y[iy]= " << y[iy] <<
"(*sp).amplitude()= " << (*sp).amplitude() << std::endl;
276 hit_signal[
chan] += ChargeFraction;
288 std::cout <<
"===========================================================" 289 "====================== " 296 std::cout <<
"end of InduceChargeFP420============================= " << std::endl;
IChargeFP420::hit_map_type induce(const CDrifterFP420::collection_type &, int numStrips, double localPitch, int numStripsW, double localPitchW, int xytype, int verbosity) override
Abs< T >::type abs(const T &t)
std::vector< AmplitudeSegmentFP420 > collection_type
static const float FiducialXYZ[3]
std::map< int, float, std::less< int > > hit_map_type
chan
lumi = TPaveText(lowX+0.38, lowY+0.061, lowX+0.45, lowY+0.161, "NDC") lumi.SetBorderSize( 0 ) lumi...
static const float capacitive_coupling[2]