StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTrsSlowAnalogSignalGenerator.cc
1 /***************************************************************************
2  *
3  * $Id: StTrsSlowAnalogSignalGenerator.cc,v 1.30 2009/11/03 14:34:19 fisyak Exp $
4  *
5  * Author:
6  ***************************************************************************
7  *
8  * Description:
9  *
10  ***************************************************************************
11  *
12  * $Log: StTrsSlowAnalogSignalGenerator.cc,v $
13  * Revision 1.30 2009/11/03 14:34:19 fisyak
14  * Remove default in zFromTB
15  *
16  * Revision 1.29 2008/10/13 19:56:12 fisyak
17  * Account that Z-offset is sector dependent
18  *
19  * Revision 1.28 2008/06/20 15:01:20 fisyak
20  * move from StTrsData to StTpcRawData
21  *
22  * Revision 1.27 2007/07/12 20:25:05 fisyak
23  * Use StarLogger, use time of flight, fix cluster shape
24  *
25  * Revision 1.26 2003/09/02 17:59:19 perev
26  * gcc 3.2 updates + WarnOff
27  *
28  * Revision 1.25 2000/06/23 00:12:41 snelling
29  * Removed dependence on local files now pointed to StDbUtilities
30  *
31  * Revision 1.24 2000/06/07 02:03:12 lasiuk
32  * exit/abort ultimatum
33  *
34  * Revision 1.23 2000/02/10 01:21:51 calderon
35  * Switch to use StTpcDb.
36  * Coordinates checked for consistency.
37  * Fixed problems with StTrsIstream & StTrsOstream.
38  *
39  * Revision 1.22 1999/12/08 02:10:42 calderon
40  * Modified to eliminate warnings on Linux.
41  *
42  * Revision 1.21 1999/10/25 18:38:49 calderon
43  * changed mPos and pos() to mPosition and position() to
44  * be compatible with StEvent/StMcEvent.
45  *
46  * Revision 1.20 1999/10/22 15:51:48 calderon
47  * Remove ifdefs for erf. Problem is solved by loading libm at the
48  * macro level.
49  *
50  * Revision 1.19 1999/10/22 00:00:15 calderon
51  * -added macro to use Erf instead of erf if we have HP and Root together.
52  * -constructor with char* for StTrsDedx so solaris doesn't complain
53  * -remove mZeros from StTrsDigitalSector. This causes several files to
54  * be modified to conform to the new data format, so only mData remains,
55  * access functions change and digitization procedure is different.
56  *
57  * Revision 1.18 1999/07/19 21:41:22 lasiuk
58  * - debug check and cleanup. No changes
59  *
60  * Revision 1.17 1999/04/27 18:33:08 lasiuk
61  * change position of time shift in sampleAnalogSignal()
62  *
63  * Revision 1.16 1999/04/27 15:05:04 lasiuk
64  * time shift in ns
65  *
66  * Revision 1.15 1999/04/23 19:19:50 lasiuk
67  * add delay to centroid of signal:
68  * Calculated in constructor (mTimeShiftOfSignalCentroid)
69  * and applied in in signalsampler()
70  *
71  * Revision 1.14 1999/02/28 20:12:50 lasiuk
72  * threshold/noise additions
73  *
74  * Revision 1.13 1999/02/26 18:26:09 lasiuk
75  * to offset must be used uniformly. Correction to
76  * "timeOfSignal" should be merged with coordinate transform
77  *
78  * Revision 1.12 1999/02/19 10:35:04 lasiuk
79  * Function prototype for StTpcCoordinateTransform
80  *
81  * Revision 1.11 1999/02/16 23:34:33 lasiuk
82  * inline 2 functions
83  * merge operations for speed up (after profiler0
84  *
85  * Revision 1.10 1999/02/15 03:38:16 lasiuk
86  * protection if min()<0
87  *
88  * Revision 1.9 1999/02/14 20:46:09 lasiuk
89  * debug info
90  *
91  * Revision 1.8 1999/02/12 01:27:18 lasiuk
92  * Limit Debug output
93  *
94  * Revision 1.7 1999/02/10 20:55:16 lasiuk
95  * Feb 10,1999
96  *
97  * Revision 1.6 1999/01/22 23:38:28 lasiuk
98  * set defaults for signal sampling and induced charge
99  *
100  * Revision 1.5 1999/01/18 21:01:30 lasiuk
101  * use fractionSampled(); enumerated types for function selection
102  *
103  * Revision 1.4 1999/01/18 10:21:10 lasiuk
104  * use integral to deposit total charge in time bin
105  *
106  * Revision 1.3 1998/11/16 14:48:19 lasiuk
107  * use wire index instead of wireNumber
108  * comment signal threshold
109  * add deltaResponse()
110  *
111  * Revision 1.2 1998/11/13 21:32:16 lasiuk
112  * adjust charge generated
113  *
114  * Revision 1.1 1998/11/10 17:12:27 fisyak
115  * Put Brian trs versin into StRoot
116  *
117  * Revision 1.5 1998/11/08 17:30:26 lasiuk
118  * allocators for SUN
119  *
120  * Revision 1.4 1998/11/04 18:47:12 lasiuk
121  * signal sampler machinery
122  *
123  * Revision 1.3 1998/10/22 14:58:27 lasiuk
124  * image charge returns double and uses PRF integral
125  *
126  * Revision 1.2 1998/10/22 00:24:27 lasiuk
127  * Oct 22
128  *
129  * Revision 1.1 1998/06/30 22:46:50 lasiuk
130  * Initial Revision
131  *
132  **************************************************************************/
133 #include <unistd.h>
134 #include <math.h>
135 #include <algorithm>
136 
137 #include "SystemOfUnits.h"
138 #include "PhysicalConstants.h"
139 #include "StDbUtilities/StTpcLocalSectorCoordinate.hh"
140 #include "StDbUtilities/StTpcPadCoordinate.hh"
141 
142 #include "StTrsSlowAnalogSignalGenerator.hh"
143 #include "TMath.h"
144 //static const double symGausAppFactor = (M_SQRT1_2*M_2_SQRTPI/(2.*mSigma1));
145 //static const double tau1 = mSigma1;
146 //static const double tau2 = 2.*mSigma1;
147 //static const double asymGausAppFactor = (M_2_SQRTPI/M_SQRT2)/(tau1+tau2);
148 //static const double asymGausUnRFactor = (M_2_SQRTPI/M_SQRT2)/(mSigma1+mSigma2);
149 // This should really come from the gas database...
150 // Hmmm...actually why should the electronics have to know about this?
151 static const double sigmaL = .05*centimeter/TMath::Sqrt(centimeter);
152 
153 StTrsAnalogSignalGenerator* StTrsSlowAnalogSignalGenerator::mInstance = 0; // static data member
154 
155 //StTrsSlowAnalogSignalGenerator::StTrsSlowAnalogSignalGenerator() {/* nopt */}
156 
157 StTrsSlowAnalogSignalGenerator::StTrsSlowAnalogSignalGenerator(StTpcGeometry* geo, StTpcSlowControl* sc, StTpcElectronics* el, StTrsSector* sec)
158  : StTrsAnalogSignalGenerator(geo, sc, el, sec)
159 {
160  // Set Defaults for the functional forms
161  mChargeDistribution = endo;
162  mSampler = symmetricGaussianApproximation;
163 
164  //
165  // Define here instead of calculating...
166  //
167 
168  mDriftVelocity = mSCDb->driftVelocity(13);
169  mSamplingFrequency = mElectronicsDb->samplingFrequency();
170 
171  mTimeBinWidth = 1./mSamplingFrequency;
172  mTau = mSigma1;
173  mTau1 = mSigma1;
174  mTau2 = 2.*mSigma2;
175  mSymGausApproxFactor = (M_SQRT1_2*M_2_SQRTPI/(2.*mSigma1));
176  mAsymGausApproxFactor = (M_2_SQRTPI/M_SQRT2)/(mTau1+mTau2);
177  mAsymGausUnRestFactor = (M_2_SQRTPI/M_SQRT2)/(mSigma1+mSigma2);
178 
179 // PR(mDriftVelocity/(centimeter/(1.e-6*second)));
180 // PR(mSamplingFrequency/MHz);
181 // PR(mTimeBinWidth/nanosecond);
182 // PR(mTau/nanosecond);
183 // PR(mTau1/nanosecond);
184 // PR(mTau2/nanosecond);
185 // PR(mSigma1/nanosecond);
186 // PR(mSigma2/nanosecond);
187 // PR(mSymGausApproxFactor);
188 // PR(mAsymGausApproxFactor);
189 // PR(mAsymGausUnRestFactor);
190 }
191 
192 StTrsSlowAnalogSignalGenerator::~StTrsSlowAnalogSignalGenerator() {/* missing */}
193 
194 StTrsAnalogSignalGenerator*
195 StTrsSlowAnalogSignalGenerator::instance()
196 {
197  if (!mInstance) {
198 #ifndef ST_NO_EXCEPTIONS
199  throw domain_error("StTrsSlowAnalogSignalGenerator::instance() Invalid Arguments");
200 #else
201  cerr << "StTrsSlowAnalogSignalGenerator::instance() Invalid Arguments" << endl;
202  cerr << "Cannot create Instance" << endl;
203  cerr << "Aborting..." << endl;
204  abort();
205 #endif
206  }
207  return mInstance;
208 }
209 
210 StTrsAnalogSignalGenerator*
211 StTrsSlowAnalogSignalGenerator::instance(StTpcGeometry* geo, StTpcSlowControl* sc, StTpcElectronics* el, StTrsSector* sec)
212 {
213  if (!mInstance) {
214  mInstance = new StTrsSlowAnalogSignalGenerator(geo, sc, el, sec);
215  }
216  // else do nothing
217  return mInstance;
218 }
219 
220 double StTrsSlowAnalogSignalGenerator::endoChargeIntegral(double xo, double yo, double xl, double xu, double yl, double yu)
221 {
222 #ifndef ST_NO_NAMESPACES
223  using namespace units;
224 #endif
225 
226  double L = (yo < mGeomDb->lastInnerSectorAnodeWire()) ?
227  mGeomDb->innerSectorAnodeWirePadPlaneSeparation() :
228  mGeomDb->outerSectorAnodeWirePadPlaneSeparation();
229  //PR(L);
230  double t3 = 1/L;
231  double t17 = -TMath::ATan(TMath::Exp(pi*(xu-xo)*t3/2))+TMath::ATan(TMath::Exp(pi*(xl-xo)*t3/2));
232  double t19 = pi*pi;
233  double t20 = 1/t19;
234  double t29 =
235  -4.0*TMath::ATan(TMath::Exp(-pi*(-yu+yo)*t3/2))*t17*t20
236  +4.0*TMath::ATan(TMath::Exp(-pi*(-yl+yo)*t3/2))*t17*t20;
237 
238  return t29;
239 }
240 
241 //double StTrsSlowAnalogSignalGenerator::gattiChargeIntegral(double xo, double yo, double xl, double xu, double yl, double yu)
242 double StTrsSlowAnalogSignalGenerator::imageChargeIntegral(double xo, double yo, double xl, double xu, double yl, double yu)
243 {
244 #ifndef ST_NO_NAMESPACES
245  using namespace units;
246 #endif
247 // cout << "StTrsSlowAnalogSignalGenerator::imageChargeIntegral" << endl;
248  // Arguments
249  // double xo, double yo, //charge location
250  // double xp, double yp, //pad centroid
251  // double xl, double xu, //integration limits
252  // double yl, double yu // "
253 
254  //
255  // Integrate a charge distribution sigma:
256  //
257  // q d
258  // sigma := 1/2 ----------------------------------
259  // 2 2 2 3/2
260  // Pi ((x - xo) + (y - yo) + d )
261  //
262  // position of charge (xo,yo)
263  // wire/pad-plane separation d
264  //
265  // yu xu
266  // / /
267  // | | q d
268  // G := | | 1/2 ---------------------------------- dx dy
269  // | | 2 2 2 3/2
270  // / / Pi ((x - xo) + (y - yo) + d )
271  // yl xl
272  //
273 
274 
275  double q = 1;
276  double d = 4*millimeter;
277 
278  //double d = (yo < mGeomDb->firstOuterSectorPadRow()) ?
279  // mGeomDb->innerSectorAnodeWirePadPlaneSeparation() :
280  // mGeomDb->outerSectorAnodeWirePadPlaneSeparation();
281  double t1 = xu*xu;
282  double t2 = xo*xo;
283  double t3 = xo*xu;
284  double t4 = t1+t2-2.0*t3;
285  double t5 = yo-yu;
286  double t8 = TMath::Power(-xu+xo,2.0);
287  double t9 = TMath::Sqrt(t8);
288  double t10 = 1/t9;
289  double t11 = yu*yu;
290  double t12 = yo*yu;
291  double t13 = yo*yo;
292  double t15 = 5*TMath::Sqrt(t1+t2+t11+(d*d)/100.-2.*t12-2.*t3+t13);
293  double t19 = TMath::ATan((50./d)*t4*t5*t10/t15);
294  double t22 = TMath::Power(-xl+xo,2.0);
295  double t23 = TMath::Sqrt(t22);
296  double t27 = xl*xl;
297  double t28 = xo*xl;
298  double t29 = t27+t2-2.0*t28;
299  double t31 = 1/t23;
300  double t33 = 5*TMath::Sqrt(t27+t2+t11+(d*d)/100-2.*t12-2.*t28+t13);
301  double t37 = TMath::ATan((50./d)*t29*t5*t31/t33);
302  double t44 = t10*t31;
303  double t46 = yo-yl;
304  double t48 = yl*yl;
305  double t49 = yo*yl;
306  double t51 = 5*TMath::Sqrt(t1+t2+t48+(d*d)/100-2.*t49-2.*t3+t13);
307  double t55 = TMath::ATan((50./d)*t4*t46*t10/t51);
308  double t62 = 5*TMath::Sqrt(t27+t2+t48+(d*d)/100-2.*t49-2.*t28+t13);
309  double t66 = TMath::ATan((50./d)*t29*t46*t31/t62);
310  double t74 = 0.1591549431*q*(-t19*xu*t23+t19*xo*t23+t37*xl*t9-t37*xo*t9)*t44
311  -0.1591549431*q*(-t55*xu*t23+t55*xo*t23+t66*xl*t9-t66*xo*t9)*t44;
312 
313  return t74;
314 }
315 
316 void StTrsSlowAnalogSignalGenerator::setChargeDistribution(StDistribution v)
317 {
318  if(v == endo ||
319  v == gatti ||
320  v == dipole )
321  mChargeDistribution = v;
322  else {
323  cerr << "Cannot Determine Charge Distribution Requested" << endl;
324  abort();
325  }
326 }
327 
328 // This is now an inline function. See the .hh file
329 // double StTrsSlowAnalogSignalGenerator::signalOnPad(double xo, double yo, double xl, double xu, double yl, double yu)
330 // {
331 
332 // // cout << "StTrsSlowAnalogSignalGenerator::signalOnPad()" << endl;
333 // switch(mChargeDistribution)
334 // {
335 // case endo:
336 // // cout << "********************Endo" << endl;
337 // return endoChargeIntegral(xo,yo,xl,xu,yl,yu);
338 // break;
339 // case gatti:
340 // // cout << "********************GATTI" << endl;
341 // // return gattiChargeIntegral(xo,yo,xl,xu,yl,yu);
342 // cout << "Gatti Distribution Not Implemented Yet!" << endl;
343 // exit(0);
344 // break;
345 // case dipole:
346 // // cout << "********************DIPOLE" << endl;
347 // return imageChargeIntegral(xo,yo,xl,xu,yl,yu);
348 // break;
349 // default:
350 // cerr << "Default Function Selected. ERROR!" << endl;
351 // exit(0);
352 // break;
353 // }
354 // }
355 
356 void StTrsSlowAnalogSignalGenerator::inducedChargeOnPad(StTrsWireHistogram* wireHistogram, Int_t sector)
357 {
358  //
359  // coordinate transform is now a data member
360  //
361  //StTpcCoordinateTransform transformer(mGeomDb, mSCDb, mElectronicsDb);
362 #if 0
363  PR(wireHistogram->minWire());
364  PR(wireHistogram->maxWire());
365 #endif
366  if(wireHistogram->minWire()<0) {
367  cerr << "Wire Plane is empty" << endl;
368  return;
369  }
370  mDriftVelocity = mSCDb->driftVelocity(sector);
371  for(int jj=wireHistogram->minWire(); jj<=wireHistogram->maxWire(); jj++) {
372 
373  // StTrsWireHistogram defines typedefs:
374  // ABOVE: typedef vector<StTrsWireBinEntry> aTpcWire
375  aTpcWire currentWire = wireHistogram->getWire(jj);
376  aTpcWire::iterator iter;
377 
378  for(iter = currentWire.begin();
379  iter != currentWire.end();
380  iter++) {
381  // What is the location of the avalanche
382  // center of Pad that is being processed?
383  // the y coordinate is the position of the wire
384 
385  //PR(*iter);
386 
387  StTpcPadCoordinate tpcRaw;
388  //StTpcLocalCoordinate xyCoord(iter->position()); // mVolumeId seems to be invalid??
389  StTpcLocalSectorCoordinate xyCoord(iter->position(),12);
390  //
391  // THIS IS IMPORTANT TO REALIZE!!!
392  //
393  // xyCoord is now:
394  // x position of hit
395  // y position of wire
396  // z drift length
397 // PR(iter->position());
398 // PR(xyCoord);
399 
400 
401 // cout << "**Transform2Raw" << endl;
402  transformer(xyCoord,tpcRaw);
403 // PR(tpcRaw);
404  int centralPad = (Int_t) tpcRaw.pad();
405  int centralRow = tpcRaw.row();
406 // PR(centralRow);
407 // PR(mDeltaRow);
408 // PR(centralPad);
409 // PR(mDeltaPad);
410 
411  // Calculate the row/pad limits
412  mRowLimits.first = (centralRow > mDeltaRow) ?
413  centralRow - mDeltaRow : centralRow;
414  mRowLimits.second = (centralRow < (mGeomDb->numberOfRows()-mDeltaRow)) ?
415  centralRow + mDeltaRow : centralRow;
416 
417  // Careful No inner/outer sector coupling!!
418  if(xyCoord.position().y() < mGeomDb->outerSectorEdge()) {
419  mRowLimits.second = min(mRowLimits.second, mGeomDb->numberOfInnerRows());
420  }
421  else {
422  mRowLimits.first = max(mRowLimits.first, (mGeomDb->numberOfInnerRows()+1));
423  mRowLimits.second = min(mRowLimits.second,(mGeomDb->numberOfRows()));
424  }
425 // PR(mRowLimits.first);
426 // PR(mRowLimits.second);
427 
428  mPadLimits.first = (centralPad > mDeltaPad) ?
429  centralPad - mDeltaPad : centralPad;
430 
431  //
432  // Loop over the cross coupled rows(irow)/pads(ipad)
433  //
434  for(int irow=mRowLimits.first; irow<=mRowLimits.second; irow++) {
435  mPadLimits.second =
436  (centralPad < (mGeomDb->numberOfPadsAtRow(irow) - mDeltaPad)) ?
437  centralPad + mDeltaPad : mGeomDb->numberOfPadsAtRow(irow);
438 // PR(mPadLimits.first);
439 // PR(mPadLimits.second);
440  for(int ipad=mPadLimits.first; ipad<=mPadLimits.second; ipad++) {
441 // cout << " row: " << irow << " pad: " << ipad << endl;
442 #ifdef ST_SECTOR_BOUNDS_CHECK
443  if( !(ipad>0 && ipad<mGeomDb->numberOfPadsAtRow(irow)) )
444  continue;
445 #endif
446  double padWidth, padLength;
447  if(irow > mGeomDb->numberOfInnerRows()) { // pad in Outer Sector
448  padWidth = mGeomDb->outerSectorPadWidth();
449  padLength = mGeomDb->outerSectorPadLength();
450  }
451  else {
452  padWidth = mGeomDb->innerSectorPadWidth();
453  padLength = mGeomDb->innerSectorPadLength();
454  }
455  tpcRaw.setPad(ipad);
456  tpcRaw.setRow(irow);
457  transformer(tpcRaw,xyCoord);
458 // PR(tpcRaw);
459 // PR(xyCoord);
460  // Integral limits for nearest pad
461  double xl = xyCoord.position().x() - padWidth/2;
462  double xu = xyCoord.position().x() + padWidth/2;
463  double yl = xyCoord.position().y() - padLength/2;
464  double yu = xyCoord.position().y() + padLength/2;
465 
466  // charge location: iter->position().x(), ycoord
467  // pad centroid: xyCoord // used to calculate integral limits
468 
469  // signalOnPad calculates charge fraction!
470  // ie-->assumes charge=1, then the total charge is scaled
471  double chargeOfSignal =
472  signalOnPad(iter->position().x(),
473  iter->position().y(), // charge location
474  xl, xu, yl, yu); // integral limits
475 // PR(chargeOfSignal);
476 // PR(iter->numberOfElectrons());
477  chargeOfSignal *= iter->numberOfElectrons();
478 // PR(chargeOfSignal);
479  //
480  // This should really be from the Coordinate transform!
481  // otherwise code has to be changed twice!
482  //
483  double timeOfSignal =
484  (iter->position().z() + mElectronicsDb->tZero()*mDriftVelocity)/mDriftVelocity;
485  // OH-OH OFFSET (replaced!...)
486  timeOfSignal =
487  iter->position().z()/mDriftVelocity;
488 
489 
490  // Check the threshold before you
491  // make and store an analog signal
492  //if() continue;
493  StTrsAnalogSignal padSignal(timeOfSignal, chargeOfSignal);
494 
495  //
496  // DIAGNOSTIC: Print out all the signals on the pad
497 // cout << "padSignal "
498 // << padSignal.time()/nanosecond << " ns\t"
499 // << padSignal.amplitude() << '\t'
500 // << irow << "," << ipad <<endl;
501  mSector->addEntry(irow,ipad,padSignal);
502  } // pad limits
503 
504  } // row limits
505 
506  } // (iterator) Loop over all wires
507 
508  } // (jj)Loop over the wires of the Histogram
509 
510 }
511 
512 /*******************************************************************/
513 // Signal Sampling
514 // For all functions below:
515 // tbin = 0,1,2,3...
516 // t = 20,150,250,350...
517 
518 double StTrsSlowAnalogSignalGenerator::deltaResponse(double tbin, StTrsAnalogSignal& s)
519 {
520  double value=0;
521 
522  // Calculate pulse Height at bin Centroid
523  double to = tbin*(1./mSamplingFrequency);
524  double deltaT = (1./mSamplingFrequency)/2;
525 
526  value = ((s.time() < to+deltaT) && (s.time() > to-deltaT)) ?
527  mGain*s.amplitude() : 0;
528 
529  value *= mFractionSampled;
530 
531  // cout << "gain/volt " << (s.amplitude()*mGain/(.001*volt)) << " mV" << endl;
532  return value;
533 }
534 
535 double StTrsSlowAnalogSignalGenerator::symmetricGaussianApproximateResponse(double tbin, StTrsAnalogSignal& s)
536 {
537  // 2
538  // 1 (t - to)
539  // F := ------------------ exp(- ---------- )
540  // 1/2 2
541  // (2 Pi) sigma 2 sigma
542 
543  // Take the value of the function at the mid-point of the
544  // time bin, and multiply it by the width of the bin!
545  // charge = F(t) dt
546  double value=0;
547 
548  // Now defined as static const double mSymGausApproxFactor
549  //double factor = 1/(TMath::Sqrt(2*pi)*mSigma1);
550  //double factor = (M_SQRT1_2*M_2_SQRTPI/(2.*mSigma1));
551 
552  // Calculate at bin Centroid (from static const double)
553  double t = mTimeBinWidth*(tbin+.5);
554 
555  value = mGain*s.amplitude()*mSymGausApproxFactor*TMath::Exp(-sqr(t - s.time())/(2*sqr(mSigma1)));
556 
557  value *= mFractionSampled*mTimeBinWidth;
558 // PR(value/(.001*volt));
559 
560  return value;
561 
562 }
563 double StTrsSlowAnalogSignalGenerator::symmetricGaussianExactResponse(double tbin, StTrsAnalogSignal& s)
564 {
565  // 2
566  // 1 (t - to)
567  // F := ------------------ exp(- ---------- )
568  // 1/2 2
569  // (2 Pi) sigma 2 sigma
570 
571  double value=0;
572 
573  // Calculate at bin Centroid
574  //double t = mTimeBinWidth*(tbin+.5);
575  double lowerBound = tbin*mTimeBinWidth; // used to be --> 100.*nanosecond;
576  double upperBound = lowerBound + mTimeBinWidth; // 100.*nanosecond;
577 // PR(lowerBound);
578 // PR(upperBound);
579 
580  value = (mGain*s.amplitude()/2.)*
581  (TMath::Erf((upperBound-s.time())/(M_SQRT2*mSigma1)) -
582  TMath::Erf((lowerBound-s.time())/(M_SQRT2*mSigma1)));
583 
584  value *= mFractionSampled;
585 
586  return value;
587 }
588 
589 double StTrsSlowAnalogSignalGenerator::asymmetricGaussianApproximateResponse(double tbin, StTrsAnalogSignal& s)
590 {
591  // Take the value of the function at the mid-point of the
592  // time bin, and multiply it by the width of the bin!
593  // charge = F(t) dt
594  double value=0;
595 
596  // Assigned as "mTau"
597  //double tau1 = mSigma1;
598  //double tau2 = 2.*mSigma1;
599 
600  // Replaced by: mAsymGausApproxFactor
601  //double factor = (M_2_SQRTPI/M_SQRT2)/(tau1+tau2);
602 
603  double t = mTimeBinWidth*(tbin+.5);
604 
605  if(t<s.time())
606  value = mGain*s.amplitude()*mAsymGausApproxFactor*TMath::Exp(-sqr(t-s.time())/(2*sqr(mTau1)));
607  else
608  value = mGain*s.amplitude()*mAsymGausApproxFactor*TMath::Exp(-sqr(t-s.time())/(2*sqr(mTau2)));
609 
610  value *= (mFractionSampled*mTimeBinWidth);
611 
612  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
613  // Amount of Charge in the bin
614  //value *= 1./mSamplingFrequency;
615 
616  // cout << "gain/volt " << (s.amplitude()*mGain/(.001*volt)) << " mV" << endl;
617  return value;
618 }
619 
620 double StTrsSlowAnalogSignalGenerator::realShaperResponse(double tbin, StTrsAnalogSignal& sig)
621 {
622  //
623  // Take the value of the function at the mid-point of the
624  // time bin, and multiply it by the width of the bin!
625  // charge = F(t) dt
626  double value=0;
627 
628  // Use mDriftVelocity...
629  //double driftVelocity = mSCDb->driftVelocity();
630 
631  // use mTau
632  //double tau = mSigma1;
633 
634  // Oh darn! Why do we need gas parmeters here...it is because
635  // we convolute the response of the electronics with the diffusion
636  // DON'T DO THAT!!!!
637  //double sigmaL = .05*centimeter/TMath::Sqrt(centimeter);
638  double t = mTimeBinWidth*(tbin+.5);
639 
640  // Remember centroid is at 2tau+10ns
641  // should be calculated via the extremem...but this
642  // can only be found numerically...There is a routine
643  // that will do this, but I won't incorporate it yet.
644  // For now use the approximation:
645  double tzero = sig.time() - 2.*mTau+10.*nanosecond;
646 
647  //double K = sigmaL*TMath::Sqrt(sig.time())/(tau*TMath::Sqrt(driftVelocity));
648  double K = sigmaL/mTau*TMath::Sqrt(tzero/mDriftVelocity);
649 
650  //double lambda = (sig.time() - t)/(K*tau) + K;
651  double lambda = (tzero - t)/(K*mTau) + K;
652 
653  // Corrected From SN197
654  value = 1./(2.*mTau)*sqr(K)*TMath::Exp(K*(lambda-.5*K))*
655  ( .5*(1+sqr(lambda))*(1-TMath::Erf(lambda/M_SQRT2)) -
656  lambda/(TMath::Sqrt(2*M_PI))*TMath::Exp(-sqr(lambda)/2));
657 
658  value *= mFractionSampled*mGain*sig.amplitude();
659 
660  // Amount of Charge in the bin
661  value *= mTimeBinWidth;
662 
663  // cout << "gain/volt " << (s.amplitude()*mGain/(.001*volt)) << " mV" << endl;
664  return value;
665 
666 }
667 
668 double StTrsSlowAnalogSignalGenerator::oneOverT(double t, double to)
669 {
670  // CAUTION: The following routine can be used to add an
671  // undershoot or unrestored baseline to the signal profile.
672  // Make sure you think you know what you are doing...
673  double value;
674 
675  if(t == to) {
676  value = 0;
677  }
678  else {
679  value = 1./(t-to); // Unrestored
680  //value = -1./(t-to); // Undershoot
681  }
682  return value;
683 }
684 
685 // double StTrsSlowAnalogSignalGenerator::asymmetricGaussianResponseWithUnrestoredBaseline(double t, StTrsAnalogSignal& sig)
686 // {
687 // double value=0;
688 
689 // // These are taken from DB at Constructor
690 // // double mSigma1 = 1;
691 // // double mSigma2 = TMath::Sqrt(8.);
692 
693 // ----> use mAsymGausUnRestFactor
694 // double factor = TMath::Sqrt(2/pi)/(mSigma1+mSigma2);
695 // if(t<sig.time())
696 // value = sig.amplitude()*mAsymGausUnRestFactor*exp(-sqr(t-sig.time())/(2*sqr(mSigma1)));
697 // else {
698 // if ((t-sig.time())<mSigma2) { // if inside 1*sigma
699 // value = sig.amplitude()*mAsymGausUnRestFactor*exp(-sqr(t-sig.time())/(2*sqr(mSigma2)));
700 // }
701 // else { // generate a fraction of undershoot/unrestored
702 // value =
703 // sig.amplitude()*mAsymGausUnRestFactor*exp(-sqr(t-sig.time())/(2*sqr(mSigma2))) +
704 // oneOverT(t,sig.time());
705 // }
706 // }
707 //
708 // value *= mFractionSampled;
709 // return (mGain*value);
710 // }
711 
712 void StTrsSlowAnalogSignalGenerator::setElectronicSampler(StSignal t)
713 {
714  if(t == delta ||
715  t == symmetricGaussianApproximation ||
716  t == symmetricGaussianExact ||
717  t == asymmetricGaussianApproximation ||
718  t == realShaper )
719  mSampler = t;
720  else {
721  cerr << "Cannot determine Electronics Response specified" << endl;
722  // this would be a good place to throw an exception
723  abort();
724  }
725 }
726 
727 // This is now an inline function. See the .hh file
728 // inline double StTrsSlowAnalogSignalGenerator::signalSampler(double t, StTrsAnalogSignal& sig)
729 // {
730 // //
731 // // This is where the function for the Signal Sampling is selected
732 // // Add a function that returns the amplitude of a signal at
733 // // a time 't' given the position in time and amplitude of all
734 // // the other signals (contained in the StTrsAnalogSignal 'sig'
735 // // -- symmetricGaussianResponse
736 // // -- asymmetricGaussianResponse
737 // // -- endoResponse
738 // if(mSampler == (StTrsSlowAnalogSignalGenerator::undefined)) {
739 // cerr << "ERROR: no function selected" << endl;
740 // exit(0);
741 // }
742 // switch(mSampler)
743 // {
744 // case symmetricGaussianApproximation:
745 // return symmetricGaussianApproximateResponse(t, sig);
746 // break;
747 // case delta:
748 // return deltaResponse(t, sig);
749 // break;
750 // case symmetricGaussianExact:
751 // return symmetricGaussianExactResponse(t, sig);
752 // break;
753 // case asymmetricGaussianApproximation:
754 // return asymmetricGaussianApproximateResponse(t, sig);
755 // break;
756 // case realShaper:
757 // return realShaperResponse(t,sig);
758 // break;
759 // //case StSignal::asymmetricGaussianResponseWithUnRestoredBaseline:
760 // //return asymmetricGaussianResponseWithUnRestoredBaseline(t, sig);
761 // //break;
762 // default:
763 // cerr << "Default Function Selected. ERROR!" << endl;
764 // exit(0);
765 // break;
766 // }
767 // }
768 
769 void StTrsSlowAnalogSignalGenerator::sampleAnalogSignal()
770 {
771  cout << "StTrsSlowAnalogSignalGenerator::sampleAnalogSignal()" << endl;
772 
773  // operates on mSector
774 
775  // I have the centroid IN TIME (make sure!!!!) of each hit!
776 #ifndef ST_NO_TEMPLATE_DEF_ARGS
777  vector<StTrsAnalogSignal> continuousAnalogTimeSequence;
778 #else
779  vector<StTrsAnalogSignal, allocator<StTrsAnalogSignal> > continuousAnalogTimeSequence;
780 #endif
781  //double freq = mElectronicsDb->samplingFrequency();
782  //PR(freq);
783  //PR(mSector->numberOfRows());
784  for(int irow=1; irow<=mSector->numberOfRows(); irow++) {
785 
786  //PR(irow);
787  for(int ipad=1; ipad<=mSector->numberOfPadsInRow(irow); ipad++) {
788  //PR(ipad);
789  continuousAnalogTimeSequence = mSector->timeBinsOfRowAndPad(irow,ipad);
790 
791  mDiscreteAnalogTimeSequence.clear();
792 
793  // Make sure it is not empty
794  //PR(continuousAnalogTimeSequence.size());
795  if(!continuousAnalogTimeSequence.size()) continue;
796  for(mTimeSequenceIterator = continuousAnalogTimeSequence.begin();
797  mTimeSequenceIterator != continuousAnalogTimeSequence.end();
798  mTimeSequenceIterator ++) {
799 // PR(mTimeSequenceIterator->time());
800 // PR(mTimeShiftOfSignalCentroid);
801  double tmpTime =
802  mTimeSequenceIterator->time() +
803  mTimeShiftOfSignalCentroid;
804  mTimeSequenceIterator->setTime(tmpTime);
805 // PR(mTimeSequenceIterator->time());
806 // PR(mTimeSequenceIterator->time()/nanosecond);
807  }
808 // cout << "row/pad " << irow << '/' << ipad << ' ' << continuousAnalogTimeSequence.size() << endl;
809 
810  // Calculate the analog signal at the centroid of the time bin
811  // Loop over all the time bins:
812 //
813 // cout << "How many signals? " << endl;
814 // PR(continuousAnalogTimeSequence.size());
815 // for(int bbb=0; bbb<continuousAnalogTimeSequence.size(); bbb++)
816 // cout << " " << bbb << " " << continuousAnalogTimeSequence[bbb] << endl;
817 // cout << "row: " << irow << " pad: " << ipad << " timeBin: " << endl;
818 
819  double timeBinT;
820  for(int itbin=0; itbin<mGeomDb->numberOfTimeBuckets(); itbin++) {
821 // cout << itbin << ", ";
822  timeBinT = itbin*mTimeBinWidth;
823 // PR(timeBinT);
824 // PR(timeBinT/nanosecond);
825 
826  double pulseHeight = 0;
827  for(mTimeSequenceIterator =continuousAnalogTimeSequence.begin();
828  mTimeSequenceIterator!=continuousAnalogTimeSequence.end();
829  mTimeSequenceIterator++) {
830 
831  //
832  // The current time bin will be filled with
833  // charge from any signal that is within
834  // 10 time bins. This should be a settable
835  // parameter.
836  //
837  if( fabs(timeBinT-mTimeSequenceIterator->time()) > 10.*mTimeBinWidth)
838  continue;
839 // cout << " tb " << itbin << " "
840 // << mTimeSequenceIterator->time()/nanosecond << " " << (*mTimeSequenceIterator) << endl;
841  pulseHeight +=
842  signalSampler(itbin, *mTimeSequenceIterator);
843  }
844 
845  //
846  // DIAGNOSTIC
847  // Print out the pulse height in each time bin
848 // cout << itbin << " pulse Height: " << pulseHeight << '\t' << (pulseHeight/(.001*volt)) << " mV" << endl;
849 
850  //
851  // Add noise here
852  //
853  // : Everywhere
854  if(!mAddNoiseUnderSignalOnly && mAddNoise) {
855  pulseHeight += generateNoise(); // noise;
856  }
857 
858  //Do not store analog Signal if it is not above a
859  // minimal threshold (should read value from database)
860 // *********************************************************
861  //if(pulseHeight < mSignalThreshold) continue;
862 // *********************************************************
863  //
864  // : Only Under Signal
865  if(mAddNoiseUnderSignalOnly && mAddNoise) {
866  //double noise = generateNoise();
867  pulseHeight += generateNoise(); //noise;
868  }
869 // cout << itbin << " pulse Height: " << pulseHeight << '\t' << (pulseHeight/(.001*volt)) << endl;
870 
871  mElectronicSignal.setTime(itbin);
872  mElectronicSignal.setAmplitude(pulseHeight);
873 // if(mElectronicSignal.amplitude() !=0 ) {
874 // if(irow == 14 && (ipad == 14 || ipad == 52)) {
875 // cout << "mElectronicSignal " << mElectronicSignal
876 // << '\t' << (mElectronicSignal.amplitude()/(.001*volt)) << endl;
877 // }
878  mDiscreteAnalogTimeSequence.push_back(mElectronicSignal);
879 
880  } // loop over time bins
881 
882  mSector->assignTimeBins(irow,ipad,mDiscreteAnalogTimeSequence);
883 
884  } // loop over pads
885 
886  } // loop over rows
887 
888 }
Generated by   doxygen 1.8.5