CMS 3D CMS Logo

InduceChargeFP420.cc
Go to the documentation of this file.
1 // File: InduceChargeFP420
3 // Date: 08.2008
4 // Description: InduceChargeFP420 for FP420
5 // Modifications:
7 //#include "SimRomanPot/SimFP420/interface/SimRPUtil.h"
9 #include <gsl/gsl_sf_erf.h>
10 #include <iostream>
11 using namespace std;
12 
13 static const float capacitive_coupling[2] = {.80, .10};
14 static const float FiducialXYZ[3] = {3.6, 4., 0.250}; // dX/2, dY/2,dZ -- mm fiducial dimension of Si plate
15 
17  int numStrips,
18  double localPitch,
19  int numStripsW,
20  double localPitchW,
21  int xytype,
22  int verbosity) {
23  signalCoupling.clear();
24  signalCoupling.push_back(capacitive_coupling[0]);
25  signalCoupling.push_back(capacitive_coupling[1]);
26 
27  // in mm local coordinates (temporarily)
28  // float FiducialX = 5., FiducialY = 10., FiducialZ = 0.250;
29  int stripLeft, stripRight, stripLeftW, stripRightW;
30  double upperBound, lowerBound, upperBoundW, lowerBoundW;
31 
32  hit_map_type hit_signal;
33 
34  // map to store pixel integrals in the x and in the y directions
35  std::map<int, float, less<int>> x, y;
36 
37  for (CDrifterFP420::collection_type::const_iterator sp = _collection_points.begin(); sp != _collection_points.end();
38  sp++) {
39  float chargePositionW = -1.; // charge in strip coord in Wide pixel
40  float chargePosition = -1.; // charge in strip coord
41 
42  // define chargePosition
43  G4ThreeVector Position3D = (*sp).position(); // charge in strip coord
44 
45  if (verbosity > 0) {
46  std::cout << " =============================*InduceChargeFP420:induce:"
47  "Position3D= "
48  << Position3D << std::endl;
49  std::cout << " localPitch= " << localPitch << std::endl;
50  std::cout << " xytype= " << xytype << " numStrips= " << numStrips << std::endl;
51  }
52 
53  // is slice still inside fiducial volume of the plate? if not -> put slice
54  // energy to zero.
55  if (abs(Position3D.x()) < FiducialXYZ[0] && abs(Position3D.y()) < FiducialXYZ[1]) {
56  if (abs(Position3D.z()) < FiducialXYZ[2]) {
57  } else {
58  // (*sp).amplitude() == 0.;
59  std::cout << " *InduceChargeFP420:Z slice outside the plate: Position3D= " << Position3D << std::endl;
60  }
61  } else {
62  // (*sp).amplitude() == 0.;
63  std::cout << " *InduceChargeFP420:XY slice outside the plate: Position3D= " << Position3D << std::endl;
64  }
65 
66  // chargePosition - still local coordinates, so exchange x and y due to 90
67  // degree rotation Yglobal::
68  if (xytype == 1) {
69  // =
70  chargePosition = 0.5 * numStrips + Position3D.x() / localPitch; // charge in strip coord. in l.r.f
71  // starting at edge of plate
72  chargePositionW = 0.5 * numStripsW + Position3D.y() / localPitchW; // charge in strip coord. in l.r.f
73  // starting at edge of plate
74  }
75  // X:
76  else if (xytype == 2) {
77  // =
78  chargePosition = 0.5 * numStrips + Position3D.y() / localPitch; // charge in strip coord. in l.r.f
79  // starting at edge of plate
80  chargePositionW = 0.5 * numStripsW + Position3D.x() / localPitchW; // charge in strip coord. in l.r.f
81  // starting at edge of plate
82  } else {
83  std::cout << "**** InduceChargeFP420: !!! ERROR: you have not to be "
84  "here !!! xytype="
85  << xytype << std::endl;
86  // break;
87  }
88 
89  if (verbosity > 0) {
90  if (xytype == 2) {
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;
96  }
97  }
98 
99  float chargeSpread = (*sp).sigma() / localPitch; // sigma in strip pitches
100  float chargeSpreadW = (*sp).sigma() / localPitchW; // sigma in strip pitches
101 
102  // Define strips intervals along x: check edge condition
103 
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);
108 
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);
113 
114  if (verbosity > 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;
122  }
125  for (int i = stripLeft; i <= stripRight; i++) {
126  /* Definition of the integration borns */
127  // go to "left"
128  if (i == 0)
129  lowerBound = 0.;
130  else {
131  gsl_sf_result result;
132  int status = gsl_sf_erf_Q_e((i - chargePosition) / chargeSpread, &result);
133  if (status != 0)
134  std::cerr << "InduceChargeFP420::could not compute gaussian tail "
135  "probability for the threshold chosen"
136  << std::endl;
137  lowerBound = 1. - result.val;
138  if (verbosity > 0) {
139  std::cout << "go to left: i= " << i << "lowerBound=" << lowerBound << std::endl;
140  }
141  }
142 
143  // go to "right"
144  if (i == numStrips - 1)
145  upperBound = 1.;
146  else {
147  gsl_sf_result result;
148  int status = gsl_sf_erf_Q_e((i - chargePosition + 1) / chargeSpread, &result);
149  if (status != 0)
150  std::cerr << "InduceChargeFP420::could not compute gaussian tail "
151  "probability for the threshold chosen"
152  << std::endl;
153  upperBound = 1. - result.val;
154 
155  if (verbosity > 0) {
156  std::cout << "go to right: i= " << i << "upperBound=" << upperBound << std::endl;
157  }
158  }
159 
160  double totalIntegrationRange = upperBound - lowerBound;
161  x[i] = totalIntegrationRange; // save strip integral
162 
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;
170  }
171 
172  } // for
175  for (int i = stripLeftW; i <= stripRightW; i++) {
176  /* Definition of the integration borns */
177  // go to "left"
178  if (i == 0)
179  lowerBoundW = 0.;
180  else {
181  gsl_sf_result result;
182  int status = gsl_sf_erf_Q_e((i - chargePositionW) / chargeSpreadW, &result);
183  if (status != 0)
184  std::cerr << "InduceChargeFP420::W could not compute gaussian tail "
185  "probability for the threshold chosen"
186  << std::endl;
187  lowerBoundW = 1. - result.val;
188 
189  if (verbosity > 0) {
190  std::cout << "go to left: i= " << i << "lowerBoundW=" << lowerBoundW << std::endl;
191  }
192  }
193 
194  // go to "right"
195  if (i == numStripsW - 1)
196  upperBoundW = 1.;
197  else {
198  gsl_sf_result result;
199  int status = gsl_sf_erf_Q_e((i - chargePositionW + 1) / chargeSpreadW, &result);
200  if (status != 0)
201  std::cerr << "InduceChargeFP420::W could not compute gaussian tail "
202  "probability for the threshold chosen"
203  << std::endl;
204  upperBoundW = 1. - result.val;
205 
206  if (verbosity > 0) {
207  std::cout << "go to right: i= " << i << "upperBoundW=" << upperBoundW << std::endl;
208  }
209  }
210 
211  double totalIntegrationRange = upperBoundW - lowerBoundW;
212  y[i] = totalIntegrationRange; // save W strip integral
213 
214  if (totalIntegrationRange <= 0.)
215  std::cout << " upperBoundW= " << upperBoundW << " lowerBoundW= " << lowerBoundW << std::endl;
216 
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;
222  }
223  }
225  // calculate signal on x strips with including capacitive coupling
226  int nSignalCoupling = signalCoupling.size();
227 
229  /*
230  int lll,copyinlll;
231  lll = unpackLayerIndex(rn0,zside);
232  if(lll==1) {copyinlll= lll/2;}
233  else if(lll==2) {copyinlll= (lll-1)/2;}
234  else{std::cout << " InduceChargeFP420:WARNING plane number in superlayer= "
235  << lll << std::endl;}
236 */
238  if (verbosity > 1) {
239  std::cout << "InduceChargeFP420: "
240  "*************************************************************** "
241  << std::endl;
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;
246  }
247  // Get the 2D charge integrals by folding x and y strips
248  for (int iy = stripLeftW; iy <= stripRightW; iy++) { // loop over Wide y index
249  for (int ix = stripLeft; ix <= stripRight; ix++) { // loop over x index
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.) {
254  // int chan = PixelDigi::pixelToChannel( ix, iy); // Get index
255  int chan = iy * numStrips + (ix + k); // Get index
256 
257  // if(k==0 ){
258  // std::cout << "InduceChargeFP420: chan= " << chan <<
259  // std::endl; std::cout << "ix= " << ix << "iy= " << iy <<
260  // std::endl;
261  // }
262  if (verbosity > 0) {
263  if (k == 0 && xytype == 2) {
264  std::cout << "InduceChargeFP420: "
265  " chan= "
266  << chan << std::endl;
267  std::cout << "ix= " << ix << "iy= " << iy << "k= " << k << "ChargeFraction= " << ChargeFraction
268  << std::endl;
269  std::cout << "hit_signal[chan]= " << hit_signal[chan] << "geVperElectron= " << geVperElectron
270  << std::endl;
271  std::cout << "signalCoupling[abs(k)]= " << signalCoupling[abs(k)] << "x[ix]= " << x[ix]
272  << "y[iy]= " << y[iy] << "(*sp).amplitude()= " << (*sp).amplitude() << std::endl;
273  }
274  }
275  // Load the amplitude:
276  hit_signal[chan] += ChargeFraction;
277  } // endif ChargeFraction
278  } // endif ix+k
279  else {
280  // std::cout << "WARNING: ix+k =" << ix+k <<
281  // std::endl;
282  } // endif ix+k
283  } // endfor k
284  } // endfor ix
285  } // endfor iy
286 
287  if (verbosity > 0) {
288  std::cout << "==========================================================="
289  "====================== "
290  << std::endl;
291  }
292 
293  } // for loop on ions (*sp)
294 
295  if (verbosity > 0) {
296  std::cout << "end of InduceChargeFP420============================= " << std::endl;
297  }
298 
299  return hit_signal;
300 }
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)
Definition: Abs.h:22
std::vector< AmplitudeSegmentFP420 > collection_type
Definition: CDrifterFP420.h:13
int k[5][pyjets_maxn]
static const float FiducialXYZ[3]
std::map< int, float, std::less< int > > hit_map_type
Definition: IChargeFP420.h:11
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]