23 int stripLeft, stripRight, stripLeftW, stripRightW;
24 double upperBound, lowerBound, upperBoundW, lowerBoundW;
30 std::map<int, float, less<int> >
x,
y;
33 for (CDrifterFP420::collection_type::const_iterator sp=_collection_points.begin(); sp != _collection_points.end(); sp++ ){
35 float chargePositionW=-1.;
36 float chargePosition=-1.;
39 G4ThreeVector Position3D = (*sp).position();
42 std::cout <<
" =============================*InduceChargeFP420:induce:Position3D= " << Position3D << std::endl;
43 std::cout <<
" localPitch= " << localPitch << std::endl;
44 std::cout <<
" xytype= " << xytype <<
" numStrips= " << numStrips << std::endl;
56 std::cout <<
" *InduceChargeFP420:Z slice outside the plate: Position3D= " << Position3D << std::endl;
61 std::cout <<
" *InduceChargeFP420:XY slice outside the plate: Position3D= " << Position3D << std::endl;
68 chargePosition = 0.5*numStrips + Position3D.x()/localPitch ;
69 chargePositionW = 0.5*numStripsW + Position3D.y()/localPitchW ;
72 else if(xytype == 2) {
74 chargePosition = 0.5*numStrips + Position3D.y()/localPitch ;
75 chargePositionW = 0.5*numStripsW + Position3D.x()/localPitchW ;
78 std::cout <<
"**** InduceChargeFP420: !!! ERROR: you have not to be here !!! xytype=" << xytype << std::endl;
84 std::cout <<
"===================================**** InduceChargeFP420: xytype= " << xytype << std::endl;
85 std::cout <<
" chargePositionW= " << chargePositionW <<
" chargePosition= " << chargePosition << std::endl;
86 std::cout <<
"Position3D= " << Position3D << std::endl;
90 float chargeSpread = (*sp).sigma()/localPitch ;
91 float chargeSpreadW = (*sp).sigma()/localPitchW ;
95 stripLeft = int( chargePosition-
clusterWidth*chargeSpread);
96 stripRight = int( chargePosition+
clusterWidth*chargeSpread);
97 stripLeft = (0<stripLeft ? stripLeft : 0);
98 stripRight = (numStrips >stripRight ? stripRight : numStrips-1);
100 stripLeftW = int( chargePositionW-
clusterWidth*chargeSpreadW);
101 stripRightW = int( chargePositionW+
clusterWidth*chargeSpreadW);
102 stripLeftW = (0<stripLeftW ? stripLeftW : 0);
103 stripRightW = (numStripsW >stripRightW ? stripRightW : numStripsW-1);
106 std::cout <<
" Position3D = " << Position3D <<
"amplitude=" << (*sp).amplitude() << std::endl;
107 std::cout <<
" MaxChargeSpread= " <<
clusterWidth*chargeSpread <<
" sigma= " << (*sp).sigma() << std::endl;
108 std::cout <<
"*** numStrips= " << numStrips <<
" localPitch= " << localPitch <<
"xytype=" << xytype << std::endl;
109 std::cout <<
" chargePosition= " << chargePosition <<
" chargeSpread= " << chargeSpread << std::endl;
110 std::cout <<
" stripLeft= " << stripLeft <<
" stripRight= " << stripRight << std::endl;
111 std::cout <<
" chargePositionW= " << chargePositionW <<
" chargeSpreadW= " << chargeSpreadW << std::endl;
112 std::cout <<
" stripLeftW= " << stripLeftW <<
" stripRightW= " << stripRightW << std::endl;
115 for (
int i=stripLeft;
i<=stripRight;
i++){
118 if (
i == 0) lowerBound = 0. ;
121 int status = gsl_sf_erf_Q_e((
i-chargePosition)/chargeSpread, &result);
122 if (status != 0)
std::cerr<<
"InduceChargeFP420::could not compute gaussian tail probability for the threshold chosen"<<std::endl;
123 lowerBound = 1. - result.val;
125 std::cout <<
"go to left: i= " <<
i <<
"lowerBound=" << lowerBound << std::endl;
130 if (
i == numStrips-1) upperBound = 1.;
133 int status = gsl_sf_erf_Q_e((
i-chargePosition+1)/chargeSpread, &result);
134 if (status != 0)
std::cerr<<
"InduceChargeFP420::could not compute gaussian tail probability for the threshold chosen"<<std::endl;
135 upperBound = 1. - result.val;
138 std::cout <<
"go to right: i= " <<
i <<
"upperBound=" << upperBound << std::endl;
142 double totalIntegrationRange = upperBound - lowerBound;
143 x[
i] = totalIntegrationRange;
145 if(totalIntegrationRange<=0.)
std::cout <<
" upperBound= " << upperBound <<
" lowerBound= " << lowerBound << std::endl;
147 std::cout <<
" *InduceChargeFP420:====================================X i = " <<
i << std::endl;
148 std::cout <<
" upperBound= " << upperBound <<
" lowerBound= " << lowerBound << std::endl;
149 std::cout <<
" totalIntegrationRange= " << totalIntegrationRange << std::endl;
150 std::cout <<
" *InduceChargeFP420:==================================== " << std::endl;
155 for (
int i=stripLeftW;
i<=stripRightW;
i++){
158 if (
i == 0) lowerBoundW = 0. ;
161 int status = gsl_sf_erf_Q_e((
i-chargePositionW)/chargeSpreadW, &result);
162 if (status != 0)
std::cerr<<
"InduceChargeFP420::W could not compute gaussian tail probability for the threshold chosen"<<std::endl;
163 lowerBoundW = 1. - result.val;
167 std::cout <<
"go to left: i= " <<
i <<
"lowerBoundW=" << lowerBoundW << std::endl;
172 if (
i == numStripsW-1) upperBoundW = 1.;
175 int status = gsl_sf_erf_Q_e((
i-chargePositionW+1)/chargeSpreadW, &result);
176 if (status != 0)
std::cerr<<
"InduceChargeFP420::W could not compute gaussian tail probability for the threshold chosen"<<std::endl;
177 upperBoundW = 1. - result.val;
180 std::cout <<
"go to right: i= " <<
i <<
"upperBoundW=" << upperBoundW << std::endl;
184 double totalIntegrationRange = upperBoundW - lowerBoundW;
185 y[
i] = totalIntegrationRange;
187 if(totalIntegrationRange<=0.)
std::cout <<
" upperBoundW= " << upperBoundW <<
" lowerBoundW= " << lowerBoundW << std::endl;
190 std::cout <<
" *InduceChargeFP420:====================================XW i= " <<
i << std::endl;
191 std::cout <<
" upperBoundW= " << upperBoundW <<
" lowerBoundW= " << lowerBoundW << std::endl;
192 std::cout <<
" totalIntegrationRange= " << totalIntegrationRange << std::endl;
193 std::cout <<
" *InduceChargeFP420:==================================== " << std::endl;
210 std::cout <<
"InduceChargeFP420: *************************************************************** " << std::endl;
211 std::cout <<
" numStripsW= " << numStripsW <<
" numStrips= " << numStrips << std::endl;
212 std::cout <<
" nSignalCoupling= " << nSignalCoupling <<
" xytype= " << xytype << std::endl;
213 std::cout <<
" stripLeftW= " << stripLeftW <<
" stripRightW= " << stripRightW << std::endl;
214 std::cout <<
" stripLeft= " << stripLeft <<
" stripRight= " << stripRight << std::endl;
217 for (
int iy=stripLeftW; iy<=stripRightW; iy++){
218 for (
int ix=stripLeft; ix<=stripRight; ix++){
219 for (
int k = -nSignalCoupling+1 ;
k<=nSignalCoupling-1 ;
k++) {
220 if (ix+
k >= 0 && ix+
k < numStrips ) {
222 if( ChargeFraction > 0. ) {
224 int chan = iy*numStrips + (ix+
k) ;
231 if(
k==0 && xytype==2){
232 std::cout <<
"InduceChargeFP420: chan= " << chan << std::endl;
233 std::cout <<
"ix= " << ix <<
"iy= " << iy <<
"k= " <<
k <<
"ChargeFraction= " << ChargeFraction << std::endl;
235 std::cout <<
"signalCoupling[abs(k)]= " <<
signalCoupling[
abs(
k)] <<
"x[ix]= " << x[ix] <<
"y[iy]= " << y[iy] <<
"(*sp).amplitude()= " << (*sp).amplitude() << std::endl;
239 hit_signal[chan] += ChargeFraction;
250 std::cout <<
"================================================================================= " << std::endl;
258 std::cout <<
"end of InduceChargeFP420============================= " << std::endl;
static float capacitive_coupling[2]
std::vector< float > signalCoupling
static float FiducialXYZ[3]
Abs< T >::type abs(const T &t)
std::map< int, float, std::less< int > > hit_map_type