StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTrsParameterizedAnalogSignalGenerator.cc
1 /***************************************************************************
2  *
3  * $Id: StTrsParameterizedAnalogSignalGenerator.cc,v 1.44 2019/02/08 21:33:54 iraklic Exp $
4  *
5  * Author: Hui Long
6  ***************************************************************************
7  *
8  * Description:tss alogrithm for signal generator in the trsRevision 1.11 2000/02/10 01:21:50 calderon
9  * normalization factors
10  ***************************************************************************
11  *
12  *
13  * $Log: StTrsParameterizedAnalogSignalGenerator.cc,v $
14  * Revision 1.44 2019/02/08 21:33:54 iraklic
15  * adding sector to the function call [iTPC related changes] : this also partial respons to ticket #3376
16  *
17  * Revision 1.43 2018/06/21 22:23:08 perev
18  * TpcGroup fixes
19  *
20  * Revision 1.42 2018/02/20 22:45:53 smirnovd
21  * Revert "Changes from Irakli's directory to make the code compile"
22  *
23  * Revision 1.40 2009/12/11 21:55:15 fisyak
24  * Account that the space between anode wires and pad plane is sensitive, bug #1715
25  *
26  * Revision 1.39 2009/11/03 14:34:19 fisyak
27  * Remove default in zFromTB
28  *
29  * Revision 1.38 2008/10/13 19:56:12 fisyak
30  * Account that Z-offset is sector dependent
31  *
32  * Revision 1.37 2008/06/20 15:01:18 fisyak
33  * move from StTrsData to StTpcRawData
34  *
35  * Revision 1.36 2007/07/12 20:25:05 fisyak
36  * Use StarLogger, use time of flight, fix cluster shape
37  *
38  * Revision 1.35 2005/12/12 21:00:12 perev
39  * 3 random generators ==> 1
40  *
41  * Revision 1.34 2005/09/09 22:12:49 perev
42  * Bug fix + IdTruth added
43  *
44  * Revision 1.32 2005/07/19 22:23:05 perev
45  * Bug fix
46  *
47  * Revision 1.31 2004/05/03 23:31:12 perev
48  * Possible non init WarnOff
49  *
50  * Revision 1.30 2004/03/26 04:18:35 fisyak
51  * Assign IdTruth even two or more tracks hit the same pixel
52  *
53  * Revision 1.29 2003/12/24 13:44:53 fisyak
54  * Add (GEANT) track Id information in Trs; propagate it via St_tpcdaq_Maker; account interface change in StTrsZeroSuppressedReaded in StMixerMaker
55  *
56  * Revision 1.28 2003/09/02 17:59:19 perev
57  * gcc 3.2 updates + WarnOff
58  *
59  * Revision 1.27 2003/05/02 23:54:22 hardtke
60  * Allow user to adjust normalFactor (i.e. Fudge Factor)
61  *
62  * Revision 1.26 2003/04/30 20:39:08 perev
63  * Warnings cleanup. Modified lines marked VP
64  *
65  * Revision 1.25 2002/04/29 21:42:15 hardtke
66  * change normalFactor back to 1.0
67  *
68  * Revision 1.24 2001/11/21 01:57:29 long
69  * adding log message for 3/2001 long;
70  * delete empty command line
71  * double InOuterFactor=1.0075,normalFactor=1.4 ;;--->
72  * double InOuterFactor=1.0075,normalFactor=1.4 ;
73  *
74  * Revision 1.23 2001/03/13 22:08:56 long
75  * *** empty log message ***
76  *
77  * Revision 1.22 2000/08/04 03:30:18 long
78  * 1)parameterize reponse function in Z
79  * 2) add normalization factor 1.4 to match real data.It should
80  * eventually go to database somewhere.
81  *
82  * Revision 1.21 2000/07/30 02:50:05 long
83  * 1) restructure code
84  * 2)add new pad response function to take into account of the charge spread
85  * in x direction additional to transverse gaussian diffusion.
86  *
87  * Revision 1.20 2000/07/17 19:19:10 long
88  * get zoffset by calling zFromTB
89  *
90  * Revision 1.19 2000/06/23 17:54:44 long
91  * offset-->offset*mTimeBinWidth
92  *
93  * Revision 1.18 2000/06/23 00:12:40 snelling
94  * Removed dependence on local files now pointed to StDbUtilities
95  *
96  * Revision 1.17 2000/06/22 17:52:44 long
97  * get zoffset by calling tBFromZ
98  *
99  * Revision 1.16 2000/06/07 02:03:11 lasiuk
100  * exit/abort ultimatum
101  *
102  * Revision 1.15 2000/05/19 17:22:29 long
103  * fix bug in central pad calculation (mCentralPad) for non_central rows.
104  *
105  * Revision 1.14 2000/04/20 21:25:20 long
106  * timeBinLowerLimit---><-----timeBinUpperLimit in
107  * "if( signalTime-timeBinT> timeBinLowerLimit*mTimeBinWidth) break;
108  * if( timeBinT-signalTime> timeBinUpperLimit*mTimeBinWidth)"
109  *
110  * Revision 1.13 2000/03/15 02:13:20 calderon
111  * Fixed bug from pad response function sigma:
112  * The data member mPadResponseFunctionSigma was assigned the right values
113  * but never used, whereas the temporary mPadRespondFunctionSigma was
114  * not initialized and then used. Removed the temporary one altogether.
115  * Also removed declaration of two_pi, use twopi from PhysicalConstants.h
116  *
117  * Revision 1.12 2000/02/24 16:35:03 long
118  * modification for step functions, normalization factors of the padresponse functions of inner ,outer sector
119  *
120  *Revision 1.12 2000/02/23 01:21:50 long
121  * modification for step functions, normalization factors of the padresponse functions of inner ,outer sector
122  * Revision 1.11 2000/02/10 01:21:50 calderon
123  * Switch to use StTpcDb.
124  * Coordinates checked for consistency.
125  * Fixed problems with StTrsIstream & StTrsOstream.
126  *
127  * Revision 1.10 2000/01/10 23:14:31 lasiuk
128  * Include MACROS for compatiblity with SUN CC5
129  *
130  * Revision 1.9 1999/11/11 19:45:11 calderon
131  * Made variables-> data members in analog signal generator to avoid
132  * initialization time when member functions are called.
133  * Inlined:
134  * StTrsParameterizedAnalogSignalGenerator::signalSampler()
135  * StTrsSector::addEntry()
136  * StTrsSector::assignTimeBins()
137  *
138  * Revision 1.8 1999/11/10 15:46:25 calderon
139  * Made changes to reduce timing, including:
140  * Made coordinate transfrom a data member of StTrsAnalogSignalGenerator
141  * Added upper-lower bound instead of symmetric cut.
142  * Revived checking if signal is above threshold.
143  *
144  * Revision 1.7 1999/11/09 19:31:33 calderon
145  * Modified loop over ContinuosAnalogTimeSequence to make it
146  * more efficient.
147  *
148  * Revision 1.6 1999/11/05 22:18:17 calderon
149  *
150  * Made private copy constructor and operator= in StTrsDigitalSector.
151  * Renamed DigitalSignalGenerators: Fast -> Old, Parameterized -> Fast
152  * and use new "Fast" as default.
153  * Added StTrsZeroSuppressedReader and StTrsZeroSuppressedReader for DAQ type
154  * data access.
155  * Removed vestigial for loop in sampleAnalogSignal() method.
156  * Write version of data format in .trs data file.
157  *
158  * Revision 1.5 1999/10/25 18:38:49 calderon
159  * changed mPos and pos() to mPosition and position() to
160  * be compatible with StEvent/StMcEvent.
161  *
162  * Revision 1.4 1999/10/22 15:51:47 calderon
163  * Remove ifdefs for erf. Problem is solved by loading libm at the
164  * macro level.
165  *
166  * Revision 1.3 1999/10/22 00:00:14 calderon
167  * -added macro to use Erf instead of erf if we have HP and Root together.
168  * -constructor with char* for StTrsDedx so solaris doesn't complain
169  * -remove mZeros from StTrsDigitalSector. This causes several files to
170  * be modified to conform to the new data format, so only mData remains,
171  * access functions change and digitization procedure is different.
172  *
173  * Revision 1.2 1999/10/06 16:50:44 long
174  * in the calculation of sigma_x,iter->position().z()----->mGeomDb->frischGrid()-iter->position().z()
175  *
176  * Revision 1.1 1999/10/04 15:43:00 long
177  * TrsParameterizedAnalogSignalGenerator using tss algorithm
178  *
179  *
180  **************************************************************************/
181 
182 #include <unistd.h>
183 #include <math.h>
184 #include <algorithm>
185 #include <ctime>
186 #include "SystemOfUnits.h"
187 #include "PhysicalConstants.h"
188 
189 #include "StTrsParameterizedAnalogSignalGenerator.hh"
190 #ifndef TPC_DATABASE_PARAMETERS
191 #include "StTpcLocalSectorCoordinate.hh"
192 #else
193 #include "StDbUtilities/StTpcLocalSectorCoordinate.hh"
194 #endif
195 #if defined (__SUNPRO_CC) && __SUNPRO_CC >= 0x500
196 using std::sort;
197 #endif
198 #include "StTrsRandom.hh"
199 #include "TMath.h"
200 static const double sigmaL = .037*centimeter/::sqrt(centimeter);
201 //static const double sigmaT = .0633*centimeter/::sqrt(centimeter);
202 static const double sqrtTwoPi = ::sqrt(twopi);
203 StTrsAnalogSignalGenerator* StTrsParameterizedAnalogSignalGenerator::mInstance = 0;
204 // static data member
205 static const double FractionCut = 0.7;// minimum Fracktion of charge to be identified as generated by the track
206 //________________________________________________________________________________
207 struct SignalSum_t {
208  double Sum;
209  int TrackId;
210 };
211 
212 //StTrsParameterizedAnalogSignalGenerator::StTrsParameterizedAnalogSignalGenerator() {/* nopt */}
213 
214 StTrsParameterizedAnalogSignalGenerator::StTrsParameterizedAnalogSignalGenerator(StTpcGeometry* geo, StTpcSlowControl* sc, StTpcElectronics* el, StTrsSector* sec)
215  : StTrsAnalogSignalGenerator(geo, sc, el, sec),
216  mPadResponseFunctionSigmaOuter(0.3913),//old 0.395
217  mPadResponseFunctionSigmaInner(0.1964)//old 0.198
218 {
219 
220  //
221  // Define here instead of calculating...
222  //
223 
224  mDriftVelocity = mSCDb->driftVelocity(13);
225 
226  mSamplingFrequency = mElectronicsDb->samplingFrequency();//hz
227  mTimeBinWidth = 1./mSamplingFrequency;//s
228  //mTau = mSigma1;
229  mTau =mElectronicsDb->tau();// s HL, 8/31/99
230  mPadResponseFunctionSigma =0; //HL,defined in member function
231  // for inner and outer setor.
232  //
233  mFractionSampled=1.0;//HL,8/31/99
234  normalFactor= 1.25; // YF from comparision AdcSum wrt data 04/04/08 | 1.0;//Default DH 5/2/03
235 
236 // PR(mDriftVelocity/(centimeter/(1.e-6*second)));
237 // PR(mSamplingFrequency/MHz);
238 // PR(mTimeBinWidth/nanosecond);
239 // PR(mTau/nanosecond);
240 
241  //
242  // Set TSS parameters for signal generation
243  //
244  mChargeFractionOuter.push_back(0.33);
245  mChargeFractionOuter.push_back(0.32);
246  mChargeFractionOuter.push_back(0.2515);
247  mChargeFractionOuter.push_back(0.0856);
248  mChargeFractionOuter.push_back(0.01607);
249  mChargeFractionOuter.push_back(0.0036);
250  mChargeFractionInner.push_back(0.333);
251  mChargeFractionInner.push_back(0.298);
252  mChargeFractionInner.push_back(0.038);
253  mChargeFractionInner.push_back(0.00181);
254  mChargeFractionInner.push_back(0.0000);
255  mChargeFractionInner.push_back(0.0000);
256 
257  mYb.push_back(0.2);
258  mYb.push_back(0.6);
259  mYb.push_back(1.0);
260  mYb.push_back(1.4);
261  mYb.push_back(1.8);
262  mNumberOfEntriesInTable=4000;
263  mRangeOfTable=4.0;
264  errorFunctionTableBuilder();
265 
266  mCentralPad = mCentralRow = 0;
267  mNumberOfRows = mGeomDb->numberOfRows();
268  mNumberOfInnerRows = mGeomDb->numberOfInnerRows();
269  mFrischGrid = mGeomDb->frischGrid();
270  rowNormalization = padWidth = padLength = 0;
271  zoffset = wire_to_plane_coupling = 0;
272  delx = gridMinusZ = sigma_x = 0;
273  dely = constant = localYDirectionCoupling = 0;
274  timeOfSignal = chargeOfSignal = 0;
275  t = tzero = K = 0;
276  sigmaLoverTau = sigmaL/mTau;
277  lambda = lambdasqr=0;
278  // srand(0);
279  localArrayBuilder();
280 
281  mAddNoise = true;
282  mNoiseRMS=1.2;
283  mAdcConversion=mElectronicsDb->adcConversion();
284  landauConstant=0.2703;
285  landauMean=2.256;
286  landauSigma=1.197;
287  expConstant=1.56538/10.;
288  expSlope=-0.589033;
289  landauCut=3.6;
290 
291  GausConstant[0]=0.2482;
292  GausMean[0]=2.317;
293  GausSigma2[0]=1.326*1.326;
294 
295  ExpConstant[0]=1.46831;
296  ExpSlope[0]=-0.569438;
297  cutT[0]=3.6;
298 
299  GausConstant[1]=0.315;
300  GausMean[1]=2.385;
301  GausSigma2[1]=1.464*1.464;
302 
303  ExpConstant[1]=1.44213;
304  ExpSlope[1]=-0.55952;
305  cutT[1]=3.7;
306 
307  GausConstant[2]=0.2153;
308  GausMean[2]=2.473;
309  GausSigma2[2]=1.636*1.636;
310 
311  ExpConstant[2]=1.4916;
312  ExpSlope[2]=-0.55763;
313  cutT[2]=4;
314 
315  GausConstant[3]=0.2;
316  GausMean[3]=2.566;
317  GausSigma2[3]=1.83*1.83;
318 
319  ExpConstant[3]=1.73945;
320  ExpSlope[3]=-0.57317;
321  cutT[3]=4.5;
322 
323  GausConstant[4]=0.1861;
324  GausMean[4]=2.599;
325  GausSigma2[4]=1.983*1.983;
326 
327  ExpConstant[4]=1.54719;
328  ExpSlope[4]=-0.544514;
329  cutT[4]=4.5;
330  GausConstant[5]=0.1734;
331  GausMean[5]=2.658;
332  GausSigma2[5]=2.17*2.17;
333 
334  ExpConstant[5]=1.70843;
335  ExpSlope[5]=-0.54852;
336  cutT[5]=5;
337 
338 
339 }
340 StTrsParameterizedAnalogSignalGenerator::~StTrsParameterizedAnalogSignalGenerator() {/* missing */}
341 
342 StTrsAnalogSignalGenerator*
343 StTrsParameterizedAnalogSignalGenerator::instance()
344 {
345  if (!mInstance) {
346 #ifndef ST_NO_EXCEPTIONS
347  throw domain_error("StTrsParameterizedAnalogSignalGenerator::instance() Invalid Arguments");
348 #else
349  cerr << "StTrsParameterizedAnalogSignalGenerator::instance() Invalid Arguments" << endl;
350  cerr << "Cannot create Instance" << endl;
351  cerr << "Exitting..." << endl;
352  exit(1);
353 #endif
354  }
355  return mInstance;
356 }
357 
358 StTrsAnalogSignalGenerator*
359 StTrsParameterizedAnalogSignalGenerator::instance(StTpcGeometry* geo, StTpcSlowControl* sc, StTpcElectronics* el, StTrsSector* sec)
360 {
361  if (!mInstance) {
362  mInstance = new StTrsParameterizedAnalogSignalGenerator(geo, sc, el, sec);
363  }
364  // else do nothing
365  return mInstance;
366 }
367 
368 void StTrsParameterizedAnalogSignalGenerator::errorFunctionTableBuilder()
369 {
370 
371  int cntr=0;
372 
373  do {
374  Double_t x = 1.0*cntr*mRangeOfTable/mNumberOfEntriesInTable;
375  Double_t y = TMath::Erf(x);
376  mErrorFunctionTable.push_back(y);
377  cntr++;
378 
379  }while(cntr < mNumberOfEntriesInTable);
380 
381 }
382 void StTrsParameterizedAnalogSignalGenerator::localArrayBuilder()
383 {
384  int rows=mNumberOfRows;
385 
386  int pad,row;
387  for(row=0;row<rows;row++)
388  {int max_pads=mGeomDb->numberOfPadsAtRow(row+1);
389  mPadsAtRow[row]=max_pads;
390  yCentroid[row]=transformer.yFromRow(20,row+1);
391  for(pad=0;pad<max_pads;pad++)
392  { xCentroid[row][pad]=transformer.xFromPad(20,row+1,pad+1);
393  // gain[row][pad]=1.0+(StTrsRandom::inst().Rndm()-0.5)*0.20;
394  }
395 
396 
397  }
398 
399 
400 }
401 
402 double StTrsParameterizedAnalogSignalGenerator::erf_fast(double argument) const
403 {
404 
405 
406  float a =argument/mRangeOfTable*mNumberOfEntriesInTable;
407  int index;
408  if(a>=0){
409  index = (int)(a+0.5);
410  if(index>=mNumberOfEntriesInTable)return 1.0;
411  return mErrorFunctionTable[index];}
412  else
413  {index = (int)(-a+0.5);
414  if(index>=mNumberOfEntriesInTable)return -1.0;
415  return -mErrorFunctionTable[index];}
416 }
417 
418 
419 
420 
421 void StTrsParameterizedAnalogSignalGenerator::inducedChargeOnPad(StTrsWireHistogram* wireHistogram, Int_t sector)
422 { double offset=transformer.zFromTB(0,sector,45);
423 
424  int PadsAtRow;
425  double sigma_xpad2;
426  double InOuterFactor=1.0075;
427  double charge_fraction[7];
428  int wire_index;
429  SignalSum_t *SignalSum = 0;//yf double *SignalSum;
430  double pitch=mGeomDb->anodeWirePitch();
431 //VPunused int PadAtRow;
432 
433  //double mPadResponseFunctionSigma;
434  //
435  // This should probably be made a data member at some point!
436 // StTpcCoordinateTransform transformer(mGeomDb, mSCDb, mElectronicsDb);
437 #if 0
438  PR(wireHistogram->minWire());
439  PR(wireHistogram->maxWire());
440 #endif
441  if(wireHistogram->minWire()<0) {
442  cerr << "Wire Plane is empty" << endl;
443  return;
444  }
445  double x,y,z,xCentroidOfPad,yCentroidOfPad;
446  double t; //VP xp,yp unused
447 //VPunused double electrons;
448  double tZero;
449  tZero=mElectronicsDb->tZero();
450  //reset all the datatbase value to keep them updated
451  mDriftVelocity = mSCDb->driftVelocity(sector);
452 
453  mSamplingFrequency = mElectronicsDb->samplingFrequency();
454  mTimeBinWidth = 1./mSamplingFrequency;
455  mTau =mElectronicsDb->tau();
456  int bin_start,bin_end;
457  int pad_start,pad_end=max(mPadsAtRow[mNumberOfInnerRows-1],mPadsAtRow[mNumberOfRows-1])+1;
458  int row_start,row_end;
459  int pad_shift,pad_shift0=-999999;
460  bin_start=0;
461  bin_end=mGeomDb->numberOfTimeBuckets();
462  row_start=0;
463  row_end= mNumberOfRows;
464 
465  int Nrbp = row_end*bin_end*pad_end;
466  SignalSum = (SignalSum_t *) calloc(Nrbp, sizeof(SignalSum_t)); assert(SignalSum);
467  //yf SignalSum=(double *)calloc(row_end*bin_end*pad_end,sizeof(double));
468  int index;
469 
470  for(int jj=wireHistogram->minWire(); jj<=wireHistogram->maxWire(); jj++) {
471 
472  // StTrsWireHistogram defines typedefs:
473  // ABOVE: typedef vector<StTrsWireBinEntry> aTpcWire
474  aTpcWire currentWire = wireHistogram->getWire(jj);
475  aTpcWire::iterator iter;
476  int Repeat;
477  // int timeBinUpperLimit = 6;
478  // int timeBinLowerLimit = 2;
479  int timeBinUpperLimit = 6;
480  int timeBinLowerLimit = 2;
481  double dAvalanch=0.015;//cm
482  int bin_low,bin_high;
483  double signalTime;
484  double pulseHeight ;
485  int max_bin=bin_end;
486  double timeBinT;
487  double *D,Dx,sigmaLz,Dt,sigmaLt;
488  double A;
489  for(iter = currentWire.begin();
490  iter != currentWire.end();
491  iter++) {
492  const int id = iter->id(); // track Id
493  z=iter->position().z();
494  D=iter->d();
495  Dx=D[0]+dAvalanch;// approximation
496  Dt=D[2]/mDriftVelocity; //mDriftVelocity : cm/s
497  // cout<<Dx<<" "<<D[2]<<" "<<mDriftVelocity<<endl;
498 
499  if( z < 0.0) { // account anode wire pad plane space
500  mCentralRow = transformer.rowFromLocal(iter->position(),sector);
501  if ((mCentralRow <= 13 && z > -0.2) ||
502  (mCentralRow > 13 && z > -0.4)) {z = - z;}
503  else {
504  cout<<"wrong z in function StTrsPara..e..."<<z<<endl;
505  continue;
506  }
507  }
508  sigmaLz=sigmaL* ::sqrt(z);
509  sigmaLt=sigmaLz/mDriftVelocity;
510  signalTime =
511  (z-offset)/mDriftVelocity; //inner outer offset is already in z
512 
513  // cout<<signalTime<<" "<<z<<endl;
514  // cin>>z;
515  if(signalTime<0.)continue;
516 
517 
518 
519  pulseHeight= iter->numberOfElectrons();
520 
521  bin_low=max(0,(int)(signalTime/mTimeBinWidth)-timeBinLowerLimit+1);
522  bin_high=min(max_bin-1,(int)(signalTime/mTimeBinWidth)+timeBinUpperLimit);
523 
524  // K = sigmaLoverTau*::sqrt(z)/mDriftVelocity;
525 
526 
527 
528  A=sigmaLt/mTau;
529  if(A<0.6){
530 
531 
532  for(int itbin=bin_low;itbin<=bin_high;itbin++){
533  SignalInTimeBin[itbin]=0.;
534  timeBinT = itbin*mTimeBinWidth;
535 
536  t=(timeBinT-signalTime)/mTau;
537  if(t<=landauCut)
538  SignalInTimeBin[itbin]=landauConstant*exp(-(t-landauMean)*(t-landauMean))/(2.*landauSigma*landauSigma)*mTimeBinWidth/mTau;
539  else
540  SignalInTimeBin[itbin]=expConstant*exp(t*expSlope)*mTimeBinWidth/mTau;
541 
542  }
543 
544  }
545  else
546  {
547  int index=(int)((A-0.6)/0.2);
548  if(index>=5.5)index=5;
549  for(int itbin=bin_low;itbin<=bin_high;itbin++){
550  SignalInTimeBin[itbin]=0.;
551  timeBinT = itbin*mTimeBinWidth;
552 
553  t=(timeBinT-signalTime)/mTau;
554  if(t<=cutT[index])
555  SignalInTimeBin[itbin]=GausConstant[index]*exp(-(t-GausMean[index])*(t-GausMean[index])/(2.*GausSigma2[index]))*mTimeBinWidth/mTau;
556  else
557  SignalInTimeBin[itbin]=ExpConstant[index]*exp(t*ExpSlope[index])*mTimeBinWidth/mTau;}
558 
559  //for
560  /* for(itbin=bin_low;itbin<=bin_high;itbin++){
561  SignalInTimeBin[itbin]=0.;
562  timeBinT = itbin*mTimeBinWidth;
563 
564  lambda = (signalTime-timeBinT)/(K*mTau) + K;
565  lambdasqr = lambda*lambda;
566  SignalInTimeBin[itbin]=0.5/mTau*(sqr(K)*exp(K*(lambda-.5*K))*
567  ( .5*(1+lambdasqr)*(1-erf_fast(lambda/M_SQRT2)) -
568 
569  lambda/sqrtTwoPi*exp(-.5*lambdasqr)))
570  *mTimeBinWidth;
571 
572  if(SignalInTimeBin[itbin]<0.)SignalInTimeBin[itbin]=0.;
573 
574  }// for itbin */
575  }//if endif
576 
577 
578 
579 
580 
581  Repeat=0;
582  x=iter->position().x();
583  y=iter->position().y();
584 
585 
586  // electrons=iter->numberOfElectrons();
587  gridMinusZ = z; // for new coordinates
588  sigma_x = iter->sigmaT();
589 
590 
591 
592 
593 
594  // StTpcLocalSectorCoordinate xyCoord(iter->position(),12);
595  // transformer(xyCoord,mTpcRaw);
596  //
597  mCentralRow = transformer.rowFromLocal(iter->position(), sector)-1;
598 
599 
600  // Calculate the row/pad limits
601  // mRowLimits.first = (mCentralRow > mDeltaRow) ?
602  //
603 
604 
605  if(y < mGeomDb->outerSectorEdge()) {
606  // mRowLimits.second = min(mRowLimits.second, mNumberOfInnerRows);
607  mRowLimits.second= mCentralRow;
608  mRowLimits.first= mCentralRow; //HL,2/20/00
609 
610  mPadResponseFunctionSigma= mPadResponseFunctionSigmaInner;
611  // rowNormalization =0.285 ; //old 1.24
612  // zoffset=mGeomDb->innerSectorzOffSet();
613  wire_to_plane_coupling=0.533*InOuterFactor;//HL,02/20/00
614  charge_fraction[0]=mChargeFractionInner[0];
615  charge_fraction[1]=mChargeFractionInner[1];
616  charge_fraction[2]=mChargeFractionInner[2];
617  charge_fraction[3]=mChargeFractionInner[3];
618  charge_fraction[4]=mChargeFractionInner[4];
619  charge_fraction[5]=mChargeFractionInner[5];
620  }
621  else {
622 
623  mRowLimits.first = max(mCentralRow - mDeltaRow, mNumberOfInnerRows);
624  mRowLimits.second = min(mCentralRow + mDeltaRow,mNumberOfRows-1);
625 
626 
627  mPadResponseFunctionSigma= mPadResponseFunctionSigmaOuter;
628  // zoffset=mGeomDb->outerSectorzOffSet();
629  // rowNormalization = 0.62;//old 2.04
630  wire_to_plane_coupling=0.512;
631  charge_fraction[0]=mChargeFractionOuter[0];
632  charge_fraction[1]=mChargeFractionOuter[1];
633  charge_fraction[2]=mChargeFractionOuter[2];
634  charge_fraction[3]=mChargeFractionOuter[3];
635  charge_fraction[4]=mChargeFractionOuter[4];
636  charge_fraction[5]=mChargeFractionOuter[5];
637  }
638 //
639 
640  //
641 
642  sigma_xpad2=sigma_x *sigma_x+ mPadResponseFunctionSigma*mPadResponseFunctionSigma;
643  for(int irow2=mRowLimits.second; irow2>=mRowLimits.first; irow2--) {
644 
645  PadsAtRow=mPadsAtRow[irow2]-1;
646  yCentroidOfPad = yCentroid[irow2];
647 
648 
649  dely = fabs(yCentroidOfPad-y);
650  wire_index=(int)(fabs(dely)/pitch+0.5);
651  if(wire_index<6)
652  localYDirectionCoupling =charge_fraction[wire_index];
653  else
654  localYDirectionCoupling =0.;
655  if(localYDirectionCoupling<1.e-5)continue;
656 
657 
658  if(Repeat<0.5){
659 
660  mCentralPad = (Int_t) transformer.padFromLocal(iter->position(),sector,irow2+1)-1;
661  if(mCentralPad> PadsAtRow)mCentralPad= PadsAtRow;//upper limit boundary check.
662 
663 
664 
665  mPadLimits.first = max(mCentralPad - mDeltaPad ,0);
666  mPadLimits.second =min(mCentralPad +mDeltaPad ,PadsAtRow);
667  pad_shift0=PadsAtRow;
668  // cout<< mPadLimits.first<< " "<<mPadLimits.second<<" "<<PadsAtRow<<" "<<irow2<<endl;
669 //
670  for(int ipad2=mPadLimits.first; ipad2<=mPadLimits.second; ipad2++) {
671 // cout << " row: " << irow << " pad: " << ipad << endl;
672 
673 
674 
675 
676 
677 
678  // Integral limits for nearest pad
679  xCentroidOfPad = xCentroid[irow2][ipad2];
680  // cout<< xCentroidOfPad << " xce "<<ipad2<<endl;
681  // cin>>hhh;
682  delx = xCentroidOfPad-x;
683 
684 
685 
686 
687  // localXDirectionCoupling[ipad2] =
688 
689  // mPadResponseFunctionSigma/::sqrt(sigma_xpad2)*exp(-0.5*delx*delx/sigma_xpad2); //::sqrt(2pi) is absorbed in the local Y coupling
690  localXDirectionCoupling[ipad2]= mPadResponseFunctionSigma/Dx*(erf_fast((Dx/2-delx)/::sqrt(2*sigma_xpad2))-erf_fast((-Dx/2-delx)/::sqrt(2*sigma_xpad2)))*0.5; //::sqrt(2pi) is absorbed in the local Y coupling
691 
692 
693  // cout<<ipad2<<" "<<localXDirectionCoupling[ipad2]<<" "<<endl;
694 
695 
696 
697 
698  // chargeOfSignal=localYDirectionCoupling*localXDirectionCoupling[ipad2]*wire_to_plane_coupling*pulseHeight*mGain*gain[irow2][ipad2];
699 
700  chargeOfSignal=localYDirectionCoupling*localXDirectionCoupling[ipad2]*pulseHeight*mGain;
701  if(chargeOfSignal<0.) continue;// ASIC threshold
702 
703 
704 
705  for(int itbin2=bin_low;itbin2<=bin_high;itbin2++){
706 
707 //VP index=irow2*bin_end*pad_end+ipad2*bin_end+itbin2;
708  index=(irow2*pad_end+ipad2)*bin_end+itbin2;
709  //yf *(SignalSum+index)+=chargeOfSignal*SignalInTimeBin[itbin2];
710  addSignal(id, chargeOfSignal*SignalInTimeBin[itbin2], *(SignalSum+index));
711 
712 
713  }
714  } // pad limits
715  Repeat=1;
716  } // if repeat<0.5
717  else{
718  pad_shift=(PadsAtRow-pad_shift0)/2;
719  // cout<<mPadLimits.first<<" "<<mPadLimits.second<<" "<<irow2<<" "<<pad_shift<<" "<<PadsAtRow<<" repeat"<<endl;
720  for(int ipad22=mPadLimits.first; ipad22<=mPadLimits.second; ipad22++) {
721 // cout << " row: " << irow << " pad: " << ipad << endl;
722 
723 
724 
725  if((ipad22+pad_shift)<0||(ipad22+pad_shift)> PadsAtRow)continue;
726 
727 
728 
729 
730 
731 
732 
733 
734 
735 
736 
737  // cout<<ipad22+pad_shift<<" repeat "<<localXDirectionCoupling[ipad22]<<" "<<endl;
738 
739  chargeOfSignal=localYDirectionCoupling*localXDirectionCoupling[ipad22]*pulseHeight*mGain ;
740  // chargeOfSignal=localYDirectionCoupling*localXDirectionCoupling[ipad22]*pulseHeight*mGain*gain[irow2][ipad22];
741 
742  if(chargeOfSignal<0.) continue;// ASIC threshold
743 
744 //
745 
746  // chargeOfSignal *= electrons*mTimeBinWidth*mGain*0.5/mTau;;
747 
748 
749 
750 
751 
752  for(int itbin22=bin_low;itbin22<=bin_high;itbin22++){
753  index=irow2*bin_end*pad_end+(ipad22+pad_shift)*bin_end+itbin22;
754 
755  //yf *(SignalSum+index)+=chargeOfSignal*SignalInTimeBin[itbin22];
756  addSignal(id, chargeOfSignal*SignalInTimeBin[itbin22], *(SignalSum+index));
757 
758  }
759 
760 
761  } // pad limits
762 
763  } //if Repeat>0.5
764  // cin>>pad_start;
765  } // row limits
766 
767 
768 
769  } // (iterator) Loop over all wires
770 
771  } // (jj)Loop over the wires of the Histogram
772 
773 
774  // cout<<" kkkkkkkkkkkk"<<endl;
775 
776  pad_start=0;
777 
778  for(int irow3=row_start;irow3<row_end;irow3++)
779  {
780 
781  int pad_end2= mPadsAtRow[irow3];
782  for(int ipad3=pad_start;ipad3<pad_end2;ipad3++)
783  { // mSector->addEntry(irow,ipad,dummy);
784  mDiscreteAnalogTimeSequence.clear();
785  for(int itbin3=bin_start;itbin3<bin_end;itbin3++)
786  {
787  index=irow3*pad_end*bin_end+ipad3*bin_end+itbin3;
788  //yf double a=*(SignalSum+index);
789  double a=(*(SignalSum+index)).Sum;
790 
791  if(a<=0)continue;
792 
793  a+=(addNoise(mNoiseRMS)*mAdcConversion);
794 
795  if(a<=mSignalThreshold)continue;
796  int id = SignalId(*(SignalSum+index));
797  mElectronicSignal.setId(id);
798  mElectronicSignal.setTime(itbin3);
799  mElectronicSignal.setAmplitude(a*normalFactor);
800  mDiscreteAnalogTimeSequence.push_back(mElectronicSignal);
801 
802  }
803 
804  if(mDiscreteAnalogTimeSequence.size()>0)mSector->assignTimeBins(irow3+1,ipad3+1,mDiscreteAnalogTimeSequence);
805  }
806  }
807  free(SignalSum);
808 
809 }
810 
811 double StTrsParameterizedAnalogSignalGenerator::realShaperResponse(double tbin, StTrsAnalogSignal& sig)
812 {
813 
814  double value=0.0;
815 
816 
817  //double sigmaL = .05*centimeter/::sqrt(centimeter);
818  // double t = mTimeBinWidth*(tbin+.5);//started from the center of time bin
819  //double t = mTimeBinWidth*(tbin);// started from the edge of time bin ,HL
820  t = tbin; //passed it directly.
821 
822  tzero = sig.time() ;// Hui Long,8/26/99
823  // tzero =0;
824  //double K = sigmaL*::sqrt(sig.time())/(tau*::sqrt(driftVelocity));
825  K = sigmaLoverTau*::sqrt((tzero- mTimeShiftOfSignalCentroid)/mDriftVelocity);//retrieve the real drift length,HL,8/31/99
826  //UNITS: sigmaL: cm/::sqrt(cm)
827  // mTau: second
828  // tzero: second
829  // mTimeBinwidth : second
830  // <mDriftVelocity :cm/second
831 
832  lambda = (tzero - t)/(K*mTau) + K;
833  lambdasqr = lambda*lambda;
834  // Corrected From SN197
835  value = .5/(mTau)*sqr(K)*exp(K*(lambda-.5*K))*
836  ( .5*(1+lambdasqr)*(1-erf_fast(lambda)) -
837 
838  lambda/sqrtTwoPi*exp(-.5*lambdasqr));
839 
840  // value*=mTimeBinWidth;
841 
842 
843  // value *= mFractionSampled*mGain*sig.amplitude();
844  value *= sig.amplitude();
845 
846  return value;
847 
848 }
849 //________________________________________________________________________________
850 void StTrsParameterizedAnalogSignalGenerator::addSignal(const int id, const double signal, SignalSum_t &S) {
851  S.Sum += signal;
852  if ( id ) {
853  if (! S.TrackId ) S.TrackId = id;
854  else // switch Id, works only for 2 tracks, more tracks ?
855  if ( S.TrackId != id && S.Sum < 2*signal) S.TrackId = id;
856  }
857  return;
858 }
859 //________________________________________________________________________________
860 int StTrsParameterizedAnalogSignalGenerator::SignalId(SignalSum_t &S) {
861  return S.TrackId > 0 ? S.TrackId : 0;
862 }
863 //________________________________________________________________________________
864 
865 double StTrsParameterizedAnalogSignalGenerator::addNoise(double sigma)
866  {
867  float x, y, z;
868 
869  y = StTrsRandom::inst().Rndm();
870  if (!y) y = StTrsRandom::inst().Rndm();
871  z = StTrsRandom::inst().Rndm();
872  x = z * 6.283185;
873 
874  return sigma*sin(x)*::sqrt(-2*::log(y));//Gaussian with mean=0
875 
876  }
877 
878 
879 void StTrsParameterizedAnalogSignalGenerator::sampleAnalogSignal()
880 {
881  cout << "StTrsParameterizedAnalogSignalGenerator::sampleAnalogSignal()" << endl;
882 
883  // operates on mSector (an StTrsSector)
884 
885  // I have the centroid IN TIME (make sure!!!!) of each hit!
886 #ifndef ST_NO_TEMPLATE_DEF_ARGS
887  vector<StTrsAnalogSignal> continuousAnalogTimeSequence;
888 //VPunused vector<StTrsAnalogSignal>::iterator lowerBound;
889 #else
890  vector<StTrsAnalogSignal, allocator<StTrsAnalogSignal> > continuousAnalogTimeSequence;
891 //VPunused vector<StTrsAnalogSignal, allocator<StTrsAnalogSignal> >::iterator lowerBound;
892 #endif
893  //double freq = mElectronicsDb->samplingFrequency();
894  //PR(freq);
895  int timeBinUpperLimit = 7;
896  int timeBinLowerLimit = 3;
897 
898  int max_pad,max_bin;
899  double timeBinT = 0;
900  int bin_low,bin_high;
901  double signalTime;
902  max_bin=mGeomDb->numberOfTimeBuckets();
903  double pulseHeight ;
904  for(int irow=1; irow<=mSector->numberOfRows(); irow++) {
905  max_pad=mSector->numberOfPadsInRow(irow);
906  for(int ipad=1; ipad<=max_pad; ipad++) {
907 
908  continuousAnalogTimeSequence = mSector->timeBinsOfRowAndPad(irow,ipad);
909 
910  mDiscreteAnalogTimeSequence.clear();
911 
912  // Make sure it is not empty
913 
914  if(!continuousAnalogTimeSequence.size()) continue;
915 
916 
917 
918 
919 
920 
921  for(mTimeSequenceIterator = continuousAnalogTimeSequence.begin();
922  mTimeSequenceIterator!=continuousAnalogTimeSequence.end();
923  mTimeSequenceIterator++) {
924 
925  signalTime = mTimeSequenceIterator->time();
926  pulseHeight=mTimeSequenceIterator->amplitude();
927  bin_low=max(0,(int)(signalTime/mTimeBinWidth)-timeBinLowerLimit+1);
928  bin_high=min(max_bin,(int)(signalTime/mTimeBinWidth)+timeBinUpperLimit);
929 
930 
931 
932  for(int itbin=bin_low;itbin<=bin_high;itbin++){
933  timeBinT = itbin*mTimeBinWidth;
934  K = sigmaLoverTau*::sqrt((signalTime- mTimeShiftOfSignalCentroid)/mDriftVelocity);
935  lambda = (-timeBinT +signalTime)/(K*mTau) + K;
936  lambdasqr = lambda*lambda;
937  SignalInTimeBin[itbin]+=(sqr(K)*exp(K*(lambda-.5*K))*
938  ( .5*(1+lambdasqr)*(1-erf_fast(lambda)) -
939 
940  lambda/sqrtTwoPi*exp(-.5*lambdasqr)))
941  *pulseHeight;
942  }// for itbin
943  } // for mtimeSequences
944 
945 
946 //
947 
948 
949 
950 
951 
952  for(int itbin=0;itbin<max_bin;itbin++){
953 //
954  // if(mAddNoise) SignalInTimeBin[itbin] += generateNoise();
955 
956  if( SignalInTimeBin[itbin]<mSignalThreshold) continue;
957 
958 
959  //Do not store analog Signal if it is not above a
960  // minimal threshold (should read value from database) rowN
961 
962 
963  mElectronicSignal.setTime(itbin);
964  mElectronicSignal.setAmplitude(SignalInTimeBin[itbin]);
965  mDiscreteAnalogTimeSequence.push_back(mElectronicSignal);
966 
967 
968  } // loop over time bins
969 //
970 
971  mSector->assignTimeBins(irow,ipad,mDiscreteAnalogTimeSequence);
972  // cout<<ipad<<" "<<irow<<endl;
973  } // loop over pads
974 
975  } // loop over rows
976 
977 
978 }
979 
980 
981 
SignalSum_t
Definition: StTpcRSMaker.h:33
Generated by   doxygen 1.8.5