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;
106 stripLeft = (0 < stripLeft ? stripLeft : 0);
107 stripRight = (numStrips > stripRight ? stripRight : numStrips - 1);
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;
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;
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;
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) {
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
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;