StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StSvtSignal.cc
1 /***************************************************************************
2  *
3  * $Id: StSvtSignal.cc,v 1.19 2015/07/29 01:48:02 smirnovd Exp $
4  *
5  * Author: Selemon Bekele
6  ***************************************************************************
7  *
8  * Description: calculates PASA signal for a gaussian input
9  *
10  ***************************************************************************
11  *
12  * $Log: StSvtSignal.cc,v $
13  * Revision 1.19 2015/07/29 01:48:02 smirnovd
14  * Added std:: to resolve ambiguity for isnan for g++ (4.8)
15  *
16  * Revision 1.18 2009/08/10 05:23:14 baumgart
17  * Adjust mPasaGain due to adjustments in initial cloud size
18  *
19  * Revision 1.17 2009/07/29 18:42:08 baumgart
20  * Increased mPasaGain to compensate for edge effects subtraction in StSvtElectronCloud.cc
21  *
22  * Revision 1.16 2009/07/02 20:29:56 baumgart
23  * Suppression of cout statements in function calcConvSignal
24  *
25  * Revision 1.15 2009/06/28 03:59:37 baumgart
26  * Compensate for addition of angular dependence in StSvtElectronCloud.cc
27  *
28  * Revision 1.14 2009/06/11 23:19:40 baumgart
29  * Small increase of mPasaGain during tune
30  *
31  * Revision 1.13 2009/03/13 22:29:37 baumgart
32  * Update mPasaGain to reflect new tune after bug fix
33  *
34  * Revision 1.12 2009/03/09 20:11:45 caines
35  * Fix to make different Rykov and Selemon methods have same gain
36  *
37  * Revision 1.11 2009/02/21 14:19:12 caines
38  * change gain to better reproduce data
39  *
40  * Revision 1.10 2008/09/22 16:03:55 caines
41  * Changed gain as needed for tuning CuCu data
42  *
43  * Revision 1.9 2005/07/23 03:37:34 perev
44  * IdTruth + Cleanup
45  *
46  * Revision 1.8 2005/02/09 14:33:35 caines
47  * New electron expansion routine
48  *
49  * Revision 1.7 2003/11/13 16:24:59 caines
50  * Further improvements to get simulator looking like reality
51  *
52  * Revision 1.6 2003/09/02 17:59:09 perev
53  * gcc 3.2 updates + WarnOff
54  *
55  * Revision 1.5 2003/07/31 19:18:10 caines
56  * Petrs improved simulation code
57  *
58  * Revision 1.4 2001/11/06 20:12:06 caines
59  * Add include for new compiler
60  *
61  * Revision 1.3 2001/08/13 15:34:18 bekele
62  * Debugging tools added
63  *
64  * Revision 1.2 2001/04/25 19:04:38 perev
65  * HPcorrs
66  *
67  * Revision 1.1 2000/11/30 20:47:49 caines
68  * First version of Slow Simulator - S. Bekele
69  *
70  **************************************************************************/
71 
72 #include "StSvtElectronCloud.hh"
73 #include "StSvtSignal.hh"
74 
75 //ClassImp(StSvtElectronCloud)
76 
77 fstream pasaOutN,pasaOut;
78 double sumAt;
79 //_______________________________________________
80 StSvtSignal::StSvtSignal()
81 {
82  mLowTBin = 1;
83  mHiTBin = 128;
84  mTotalHitCharge = 0.0;
85  mFractionOfCharge = 0.0;
86  mCollectedCharge = 0.0;
87 
88  mSigmaMajor = 0.0;
89  mSigmaMinor = 0.0;
90  mDriftVel = 0.0; //[in mm/micro seconds];
91  mTimeCenter = 0.0;
92  mTimeWidth = 0.0;
93 
94  GAP_TWIDTH = 36*1.e-6;
95 
96  memset(mSignal,0,sizeof(mSignal[0])*128);
97 
98  // mPasaGain = 7.2; // uV/e Original number
99  mPasaGain = 13.4; // Continuing tuning
100 
101  for(int i = 0; i < 4; i++)
102  mPasa[i] = 0.0;
103 
104  mTau_s = 11.0*0.001; //in micro seconds
105  mTau_l = 500.0*0.001; //in micro seconds
106 
107 }
108 //_______________________________________________
109 StSvtSignal::~StSvtSignal()
110 {
111 
112 }
113 //_______________________________________________
114 void StSvtSignal::setAnodeTimeBinSizes(double timeBinSize,double anodeSize)
115 {
116  mTimeBinSize = timeBinSize;
117  mAnodeSize = anodeSize;
118 
119 }
120 
121 //_______________________________________________
122 void StSvtSignal::setDriftVelocity(double driftVelocity)
123 {
124 mDriftVel = driftVelocity;
125 }
126 
127 //_______________________________________________
128 void StSvtSignal::setOption(int option)
129 {
130  mOption = option;
131 }
132 
133 //_______________________________________________
134 void StSvtSignal::pasaRelatedStuff()
135 {
136 
137  //**** pasa related stuff
138 
139  //for Rykov's piece of code
140 
141  unNormPasaConst();
142 
143  //To speed things up do not calculate these values every time, just use the already calculated values
144 
145  mPeakTimeR = 0.04308;
146  mPasaMaxR = 4.4026;
147  mFwhmR = 0.049364;
148 
149  /*
150  peakingTimeR(); // [micro seconds]
151  cout<<"mPeakTimeR = "<<mPeakTimeR<<endl;
152  cout<<"mPasaMaxR = "<<mPasaMaxR<<endl;
153  halfWidthAtHalfMaxR();
154  cout<<"mFwhmR = "<<mFwhmR<<endl;
155  */
156 
157  normPasaConst();
158  arrays();
159 
160  //for Selemons piece of code
161  //To speed things up do not calculate these values every time, just use the already calculated values
162 
163  mPeakTimeS = 0.04308;
164  mPasaMaxS = 2.74619e-09;
165  mFwhmS = 0.0494;
166  // mPasaNorm = 2.62181e+09; // [micro volts]/[micro seconds**4]-e
167  mPasaNorm = mPasaGain*3.64140278e+08; // Automatic Conversion
168 
169  /*
170  peakingTimeS();
171  cout<<"\nmPeakTimeS = "<<mPeakTimeS<<endl;
172  cout<<"mPasaMaxS = "<<mPasaMaxS<<endl;
173  halfWidthAtHalfMaxS();
174  cout<<"mFwhmS = "<<mFwhmS<<endl;
175  mPasaNorm = mPasaGain/mPasaMaxS; // [micro volts]/[micro seconds**4]-e
176  cout<<"mPasaNorm = "<<mPasaNorm<<endl;
177  */
178 }
179 
180 //_______________________________________________
181 void StSvtSignal::doPasaOnly(int option){
182  double t = 0.0;
183  // double tStep = 0.04;
184  double tStep = mTimeBinSize*0.001; //microseconds
185  double pasaFunValue;
186 
187  if(option)
188  pasaOutN.open("pasaNorm.dat",ios::out);
189  else
190  pasaOut.open("pasa.dat",ios::out);
191 
192  unNormPasaConst();
193 
194  for(int j = 0; j <25000; j++){
195 
196  t = t + tStep;
197 
198  pasaFunValue = pasaRes(t);
199 
200  if(option){
201  pasaFunValue = mPasaNorm*pasaFunValue;
202  pasaOutN<<t<<setw(20)<<pasaFunValue<<endl;
203 
204  } else {
205 
206  pasaOut<<t<<setw(20)<<pasaFunValue<<endl;
207  }
208  /*
209  if(j < 100){
210  cout<<t<<"\t"<<pasaFunValue<<endl;
211  }
212  */
213  }
214 
215  halfWidthAtHalfMaxS();
216 
217  if(option)
218  pasaOutN.close();
219  else
220  pasaOut.close();
221 
222 }
223 //_______________________________________________
224 void StSvtSignal::setCloud(StSvtElectronCloud* elCloud)
225 {
226  double sigmaSqDiff,sigmaSqDiff2;
227  double cosPhi,cosPhi2,sinPhi,sinPhi2;
228 
229  /* old stuff
230  sigmaMajor2 = elCloud->getSigma1Sq();
231  sigmaMinor2 = elCloud->getSigma2Sq();
232 
233  sigmaMajor = ::sqrt(sigmaMajor2); // [mm]
234  sigmaMinor = ::sqrt(sigmaMinor2); // [mm]
235  sigmaSqDiff = sigmaMajor2 - sigmaMinor2;
236  sigmaSqDiff2 = sigmaSqDiff*sigmaSqDiff;
237 
238  mSigmaMajor = elCloud->getSigma1(); // [mm]
239  mSigmaMinor = elCloud->getSigma2(); // [mm]
240  //mSigmaMinor = (mSigmaMinor/0.24)*(4.78/128); // [micro seconds]
241  //mSigmaMinor = mSigmaMinor*0.0001555989583; // [micro seconds]
242  */
243 
244  mTrackId = elCloud->getTrackId();
245  mSigmaMajor = elCloud->getSigmaDrift();
246  double sigmaMajor2 = mSigmaMajor*mSigmaMajor;
247  mSigmaMinor = elCloud->getSigmaAnode();
248  double sigmaMinor2 = mSigmaMinor*mSigmaMinor;
249  sigmaSqDiff = sigmaMajor2 - sigmaMinor2;
250  sigmaSqDiff2 = sigmaSqDiff*sigmaSqDiff;
251 
252  mTotalHitCharge = elCloud->getChargeAtAnode();
253 
254  cosPhi = cos(elCloud->getPhi());
255  sinPhi = sin(elCloud->getPhi());
256  cosPhi2 = cosPhi*cosPhi;
257  sinPhi2 = sinPhi*sinPhi;
258 
259  mC1 = (-0.39894228*cosPhi*sinPhi*sigmaSqDiff)/mSigmaMajor; // [mm]
260  mC2 = (0.39894228*cosPhi2*sinPhi2*sigmaSqDiff2)/sigmaMajor2; // [mm]**2
261  mC3 = (sigmaMajor2*sigmaMinor2/sigmaMajor2) + (cosPhi2*sinPhi2*sigmaSqDiff2/sigmaMajor2); // [mm]**2
262 }
263 
264 //_______________________________________________
265 double StSvtSignal::chargeFraction(int an, double anHit)
266 {
267  //anodes are numbered from 1
268  double anodeHitCent = anHit*mAnodeSize;
269 
270  mAnRightEdge = an*mAnodeSize;
271  mAnLeftEdge = (an - 1)*mAnodeSize;
272 
273  //Fraction of total charge collected
274 
275  mFractionOfCharge = prob1(mAnRightEdge - anodeHitCent,mSigmaMajor) - prob1(mAnLeftEdge - anodeHitCent,mSigmaMajor);
276 
277  //cout<<"mAnRightEdge = "<<mAnRightEdge<<endl;
278  //cout<<"mAnLeftEdge = "<<mAnLeftEdge<<endl;
279  //cout<<"anodeHitCent = "<<anodeHitCent<<endl;
280  //cout<<"mSigmaMajor = "<<mSigmaMajor<<endl;
281 
282  if(mFractionOfCharge < 0.000001)
283  {
284  mCollectedCharge = 0.0;
285  }
286  else
287  {
288  //cout<<"mFractionOfCharge = "<<mFractionOfCharge<<endl;
289  mCollectedCharge = mTotalHitCharge*mFractionOfCharge;
290  }
291 
292  return mCollectedCharge;
293 
294 }
295 
296 //_______________________________________________
297 //relTimeCenter = Center of gravity relative to the cloud y-center of gravity
298 
299 int StSvtSignal::timeCenterAndWidth(double anHit,double timeHit)
300 {
301  double leftArgument,leftArgument2,rightArgument,rightArgument2;
302  double exp1,exp2,relTimeCenter,relTimeCenter2,timeWidth2,driftVel2;
303 
304  double anodeHitCent = anHit*mAnodeSize;
305  double mTimeHit = timeHit*mTimeBinSize;
306 
307  leftArgument = (mAnLeftEdge - anodeHitCent)/mSigmaMajor;
308  leftArgument2 = leftArgument*leftArgument;
309  rightArgument = (mAnRightEdge - anodeHitCent)/mSigmaMajor;
310  rightArgument2 = rightArgument*rightArgument;
311 
312  //cout<<"leftArgument = "<<leftArgument<<endl;
313  //cout<<"rightArgument = "<<rightArgument<<endl;
314 
315  exp1 = exp(-0.5*leftArgument2);
316  exp2 = exp(-0.5*rightArgument2);
317 
318  //if(exp1 < 1e-30)
319  // exp1 = 0.0;
320  //if(exp2 < 1e-30)
321  //exp2 = 0.0;
322 
323  //cout<<"exp1 = "<<exp1<<endl;
324  //cout<<"exp2 = "<<exp2<<endl;
325  //cout<<"mC1 = "<<mC1<<endl;
326  //cout<<"mC2 = "<<mC2<<endl;
327  //cout<<"mC3 = "<<mC3<<endl;
328 
329  relTimeCenter = mC1*(exp1 - exp2); //mm
330  //cout<<"relTimeCenter = "<<relTimeCenter<<endl;
331  timeWidth2 = mC2*(leftArgument*exp1 - rightArgument*exp2) + mC3*mFractionOfCharge; //[mm]**2
332 
333  // center of gravity
334  relTimeCenter = relTimeCenter/mFractionOfCharge;
335  //cout<<"relTimeCenter = "<<relTimeCenter<<endl;
336  relTimeCenter2 = relTimeCenter*relTimeCenter;
337  timeWidth2 = timeWidth2/mFractionOfCharge;
338  timeWidth2 = fabs(timeWidth2 - relTimeCenter2);
339  //cout<<"timeWidth2 = "<<timeWidth2<<endl;
340 
341  driftVel2 = mDriftVel*mDriftVel;
342 
343  mTimeCenter = mTimeHit + (relTimeCenter/mDriftVel); //micro seconds
344  //cout<<" mTimeCenter = "<< mTimeCenter<<endl;
345  mTimeWidth = ::sqrt((timeWidth2/driftVel2) + GAP_TWIDTH); //micro seconds
346  //cout<<"mTimeWidth = "<<mTimeWidth<<endl;
347 
348  if(std::isnan(mTimeWidth))
349  return 1;
350  else return 0;
351 
352 }
353 
354 //_______________________________________________
355 void StSvtSignal::resetPeakAndUnderShoot()
356 {
357  mPeakSignal = 0.0;
358  mMinUnderShoot = 0.0;
359 }
360 //_______________________________________________
361 void StSvtSignal::setTimeWidth(double timWidth)
362 {
363  mTimeWidth = timWidth;
364  //cout<<"mTimeWidth = "<<mTimeWidth<<endl;
365 }
366 
367 //_______________________________________________
368 void StSvtSignal::calcConvSignal(double chargeOnAnode)
369 {
370  int nMin, nMax;
371  double tStep = 0;
372 
373  resetSignal(mLowTBin,mHiTBin);
374  resetPeakAndUnderShoot();
375 
376  //4.78/128 = 0.03734375;
377  tStep = mTimeBinSize; //microseconds
378 
379  //nMin = (int)((mTimeCenter - 4*mTimeWidth)/tStep);
380  int mTCenter = (int)(mTimeCenter/tStep);
381  nMin = mTCenter - 10;
382  if(nMin <= 0) nMin = 1;
383 
384  mLowTBin = nMin;
385 
386  //nMax = (int)((mTimeCenter + 5*mTimeWidth)/tStep);
387  nMax = mTCenter + 20;
388  if(nMax > 128) nMax = 128;
389 
390  mHiTBin = nMax;
391 
392  if(mOption == 1)
393  {
394  //cout<<"using Rykove's method"<<endl;
395  rykovSignal(nMin,nMax, tStep);
396  }
397  else if(mOption == 2 )
398  {
399  //cout<<"using Selemon's version"<<endl;
400  selemonSignal(nMin,nMax,tStep,chargeOnAnode);
401  }
402  else
403  {
404  //cout<<"using both versions"<<endl;
405  if(mTimeWidth< 0.14)
406  //if(mTimeWidth > 0.02 && mTimeWidth< 0.14)
407  {
408  //cout<<"now using Rykove's version"<<endl;
409  //cout << "Rykov" << mTCenter << endl;
410  rykovSignal(nMin,nMax, tStep);
411  }
412  else
413  {
414  //cout<<"now using selemons version"<<endl;
415  //cout << "Selemon" << mTCenter << endl;
416  selemonSignal(nMin,nMax,tStep,chargeOnAnode);
417  }
418  }
419 
420 }
421 //_______________________________________________
422 void StSvtSignal::rykovSignal(int nMin,int nMax, double tStep)
423 {
424  double t = 0;
425 
426  for(int n = nMin; n <= nMax ; n++)
427  {
428  t = n*tStep;
429  mSignal[n - 1] = signal(t)*1000000.0; // [micro volts]
430  if(n < (int)(mTimeCenter/tStep) && mSignal[n - 1] < 0.0)
431  mSignal[n - 1] = 0;
432  if(mSignal[n - 1] > mPeakSignal)
433  mPeakSignal = mSignal[n - 1];
434  if(mSignal[n - 1] < mMinUnderShoot)
435  mMinUnderShoot = mSignal[n - 1];
436  }
437 }
438 //_______________________________________________
439 void StSvtSignal::selemonSignal(int nMin,int nMax, double tStep, double charge)
440 {
441  double t = 0, sig = 0;
442  int numOfIntPoints = 0;
443  sumAt = 0.0;
444 
445  for(int n = nMin; n <= nMax ; n++)
446  {
447  t = n*tStep;
448  if(mTimeWidth <= 0.06)
449  {
450  sig = analConvInt(t,mTimeWidth,mTimeCenter);
451  mSignal[n - 1] = charge*mPasaNorm*sig; // [micro volts]
452  if(n < (int)(mTimeCenter/tStep) && mSignal[n - 1] < 0.0)
453  mSignal[n - 1] = 0;
454  if(mSignal[n - 1] > mPeakSignal)
455  mPeakSignal = mSignal[n - 1];
456  if(mSignal[n - 1] < mMinUnderShoot)
457  mMinUnderShoot = mSignal[n - 1];
458  }
459  else
460  {
461  numOfIntPoints = 2;
462  sig = numConvInt(nMin, n, numOfIntPoints, tStep, t);
463  //cout<<"sig = "<<sig<<endl;
464  mSignal[n - 1] = charge*mPasaNorm*sig; // [micro volts]
465  if(n < (int)(mTimeCenter/tStep) && mSignal[n - 1] < 0.0)
466  mSignal[n - 1] = 0;
467  if(mSignal[n - 1] > mPeakSignal)
468  mPeakSignal = mSignal[n - 1];
469  if(mSignal[n - 1] < mMinUnderShoot)
470  mMinUnderShoot = mSignal[n - 1];
471  }
472  }
473 
474 }
475 
476 //_______________________________________________
477 //pasa dependent coefficients
478 void StSvtSignal::unNormPasaConst()
479 {
480 
481  double p = mTau_s/mTau_l;
482  double q = 1.0/(1.0 - p);
483  mPasa[4] = p;
484 
485  for(int i = 4; i > 0; i--)
486  {
487  mPasa[i-1] = i*mPasa[i]*q;
488  }
489  mPasa[4] = 1.0;
490 
491 }
492 //_______________________________________________
493 void StSvtSignal::peakingTimeS()
494 {
495 
496  double t = 0.0;
497  // double tStep = 0.04;
498  double tStep = mTimeBinSize*0.001; //microseconds
499  double pasaFunValue = -1.e20;
500  mPasaMaxS = -2.0e20;
501 
502  do {
503  mPasaMaxS = pasaFunValue;
504  t = t + tStep;
505 
506  pasaFunValue = pasaRes(t);
507  //cout<<"t = "<<t<<"\tpasaFunValue = "<<pasaFunValue<<endl;
508 
509  } while( pasaFunValue > mPasaMaxS);
510 
511  mPeakTimeS = t - tStep;
512 
513 }
514 //_______________________________________________
515 void StSvtSignal::peakingTimeR()
516 {
517 
518  double t = 0.0;
519  // double tStep = 0.04;
520  double tStep = (mTimeBinSize*0.001)/mTau_s; //microseconds
521  double pasaFunValue = -1.e20;
522  mPasaMaxR = -2.0e20;
523 
524  do {
525  mPasaMaxR = pasaFunValue;
526  t = t + tStep;
527  pasaFunValue = mPasa[0];
528  double c1 = 1.0;
529 
530  for(int i = 1; i <= 4; i++)
531  {
532  c1 = c1*t;
533  pasaFunValue = pasaFunValue + mPasa[i]*c1;
534  }
535 
536  pasaFunValue = pasaFunValue*exp(-t) - mPasa[0]*exp(-(mTau_s/mTau_l)*t);
537  //cout<<"t = "<<t<<"\tpasaFunValue = "<<pasaFunValue<<endl;
538 
539  } while( pasaFunValue > mPasaMaxR);
540 
541  double peakTime = t - tStep;
542  mPeakTimeR = mTau_s*peakTime;
543 
544 }
545 //_______________________________________________
546 void StSvtSignal::halfWidthAtHalfMaxS()
547 {
548  double t;
549  //double tStep = 0.04;
550  double tStep = mTimeBinSize*0.001; //microseconds
551  double pasaFunValue = 0.0;
552  mFwhmS = 0.0;
553 
554  for(int i = 0; i < 2; i++)
555  {
556  mFwhmS = -1.0*mFwhmS;
557  t = mPeakTimeS;
558  pasaFunValue = mPasaMaxS; // from previous function
559 
560  do {
561  t = t - tStep;
562  pasaFunValue = pasaRes(t);
563 
564  } while( pasaFunValue > 0.5*mPasaMaxS);
565 
566  mFwhmS = mFwhmS + t;
567  tStep = -tStep;
568  }
569 
570 }
571 //_______________________________________________
572 void StSvtSignal::halfWidthAtHalfMaxR()
573 {
574  double t;
575  //double tStep = 0.04;
576  double tStep = mTimeBinSize*0.001; //microseconds
577  double pasaFunValue = 0.0;
578  mFwhmR = 0.0;
579 
580  for(int i = 0; i < 2; i++)
581  {
582  mFwhmR = -1.0*mFwhmR;
583  t = mPeakTimeR/mTau_s;
584  pasaFunValue = mPasaMaxR; // from previous function
585 
586  do {
587  t = t - tStep;
588  pasaFunValue = mPasa[0];
589  double c1 = 1.0;
590 
591  for(int i = 1; i <= 4; i++)
592  {
593  c1 = c1*t;
594  pasaFunValue = pasaFunValue + mPasa[i]*c1;
595  }
596 
597  pasaFunValue = pasaFunValue*exp(-t) - mPasa[0]*exp(-(mTau_s/mTau_l)*t);
598 
599  } while( pasaFunValue > 0.5*mPasaMaxR);
600 
601  mFwhmR = mFwhmR + t;
602  tStep = -tStep;
603  }
604 
605  mFwhmR = mTau_s*mFwhmR;
606 }
607 //_______________________________________________
608 void StSvtSignal::normPasaConst()
609 {
610  double s1 = mPasaGain*1.e-6/mPasaMaxR;
611 
612  for(int i = 0; i <= 4; i++)
613  mPasa[i] = s1*mPasa[i];
614 
615  //p1_ras1 = mTau_l*mPasa[0];
616 
617 }
618 
619 
620 void StSvtSignal::arrays()
621 {
622  mArray1[0] = 2.46196981473530512524e-10;
623  mArray1[1] = 0.564189564831068821977;
624  mArray1[2] = 7.46321056442269912687;
625  mArray1[3] = 48.6371970985681366614;
626  mArray1[4] = 196.520832956077098242;
627  mArray1[5] = 526.445194995477358631;
628  mArray1[6] = 934.528527171957607540;
629  mArray1[7] = 1027.55188689515710272;
630  mArray1[8] = 557.535335369399327526;
631 
632  mArray2[0] = 1.0;
633  mArray2[1] = 13.2281951154744992508;
634  mArray2[2] = 86.7072140885989742329;
635  mArray2[3] = 354.937778887819891062;
636  mArray2[4] = 975.708501743205489753;
637  mArray2[5] = 1823.90916687909736289;
638  mArray2[6] = 2246.33760818710981792;
639  mArray2[7] = 1656.66309194161350182;
640  mArray2[8] = 557.535340817727675546;
641 
642  mArray3[0] = 0.0;
643  mArray3[1] = 0.564189583547755073984;
644  mArray3[2] = 1.275366707599781044160;
645  mArray3[3] = 5.019050422511804774140;
646  mArray3[4] = 6.160210979930535851950;
647  mArray3[5] = 7.409742699504489391600;
648  mArray3[6] = 2.978866653721002406700;
649 
650  mArray4[0] = 1.0;
651  mArray4[1] = 2.26052863220117276590;
652  mArray4[2] = 9.39603524938001434673;
653  mArray4[3] = 12.0489539808096656605;
654  mArray4[4] = 17.0814450747565897222;
655  mArray4[5] = 9.60896809063285878198;
656  mArray4[6] = 3.36907645100081516050;
657 
658  mArray5[0] = 0.0;
659  mArray5[1] = 9.60497373987051638749;
660  mArray5[2] = 90.0260197203842689217;
661  mArray5[3] = 2232.00534594684319226;
662  mArray5[4] = 7003.32514112805075473;
663  mArray5[5] = 55592.3013010394962768;
664 
665  mArray6[0] = 1.0;
666  mArray6[1] = 33.5617141647503099647;
667  mArray6[2] = 521.357949780152679795;
668  mArray6[3] = 4594.32382970980127987;
669  mArray6[4] = 22629.0000613890934246;
670  mArray6[5] = 49267.3942608635921086;
671 
672 }
673 //_______________________________________________
674 double StSvtSignal::signal(double t)
675 {
676  double signal;
677  double localTime = 0.0;
678 
679  localTime = t - mTimeCenter;
680 
681  if(mCollectedCharge < 1.0)
682  signal = 0;
683  else if(localTime + 5.0*mTimeWidth < 0.0)
684  signal = 0;
685  else if(mTimeWidth <= 0.02*mTau_s)
686  signal = getShortSignal(localTime);
687  else
688  signal = getLongSignal(localTime);
689 
690  return signal;
691 
692 }
693  //_______________________________________________
694 double StSvtSignal::getShortSignal(double localTime)
695  {
696  double signal, sig_s, sig_l;
697 
698  signal = -mPasa[0]*exp(-localTime/mTau_l);
699  localTime = localTime/mTau_s;
700  sig_s = 0.0;
701  sig_l = 1.0;
702 
703  for(int i = 0; i <= 4; i++)
704  {
705  sig_s += mPasa[i]*sig_l;
706  sig_l = sig_l*localTime;
707  }
708 
709  signal = mCollectedCharge*( signal + sig_s*exp(-localTime));
710 
711  return signal;
712  }
713 
714 //_______________________________________________
715 double StSvtSignal::getLongSignal(double localTime)
716  {
717  double ds0, ds1, ds2, dsc, dsc0, dsc1, dsc2, dsc3;
718  double da1,da2,a2,a3,a4;
719 
720  da1 = mTimeWidth/mTau_s;
721  da2 = da1*da1;
722  a2 = 0.5*mCollectedCharge/da1;
723  a3 = mTimeWidth/mTau_l;
724  a4 = 0.5*a3*a3;
725 
726  ds0 = localTime/mTimeWidth;
727  ds1 = ds0 - da1;
728  ds2 = da1*ds1;
729  dsc = fabs(0.7071067811865476*ds1);
730  dsc2 = dsc*dsc;
731  dsc3 = exp(-0.5*ds0*ds0);
732 
733  if(dsc < 1.0)
734  dsc0 = useArrays5And6(ds1,dsc);
735  else
736  dsc0 = useArrays3And4Or1And2(ds1,dsc);
737 
738  dsc1 = 0.3989422805274159*da1;
739  if(dsc >= 1.0 && ds1 < 0.0)
740  dsc1 = dsc1 + ds2*dsc0;
741  else
742  {
743  dsc0 = dsc0*exp(-da1*(ds1 + 0.5*da1));
744  dsc1 = dsc1*dsc3 + ds2*dsc0;
745  }
746 
747  double sig0 = mPasa[0]*dsc0 + mPasa[1]*dsc1;
748  for(int i = 2; i <= 4; i++)
749  {
750  dsc2 = da2*(i - 1)*dsc0 + ds2*dsc1;
751  sig0 = sig0 + mPasa[i]*dsc2;
752  dsc0 = dsc1;
753  dsc1 = dsc2;
754  }
755 
756  if(dsc >= 1.0 && ds1 < 0.0)
757  sig0 = sig0*dsc3;
758 
759  double sig1 = ds0 - a3;
760 
761  sig0 = sig0 - mPasa[0]*exp(a4 - localTime/mTau_l)*prob1(sig1,1.0);
762  double signal = mCollectedCharge*sig0;
763 
764  return signal;
765 }
766 //_______________________________________________
767 double StSvtSignal::useArrays5And6(double ds1, double dsc)
768  {
769  double sigUp,sigDn;
770 
771  sigUp = mArray5[0];
772  sigDn = mArray6[0];
773 
774  for(int i = 1; i <= 5; i++)
775  {
776  sigUp = sigUp*dsc*dsc + mArray5[i];
777  sigDn = sigDn*dsc*dsc + mArray6[i];
778  }
779 
780  double dsc0 = 0.5 - 0.5*dsc*sigUp/sigDn;
781  if(ds1 > 0.0)
782  dsc0 = 1.0 - dsc0;
783 
784  return dsc0;
785 }
786 //_______________________________________________
787 double StSvtSignal::useArrays3And4Or1And2(double ds1,double dsc)
788  {
789  double sigUp,sigDn;
790 
791  if(dsc >= 8.0)
792  {
793  sigUp = mArray3[0];
794  sigDn = mArray4[0];
795 
796  for(int i = 1; i <= 6; i++)
797  {
798  sigUp = sigUp*dsc + mArray3[i];
799  sigDn = sigDn*dsc + mArray4[i];
800  }
801  }
802  else
803  {
804  sigUp = mArray1[0];
805  sigDn = mArray2[0];
806 
807  for(int i = 1; i <= 8; i++)
808  {
809  sigUp = sigUp*dsc + mArray1[i];
810  sigDn = sigDn*dsc + mArray2[i];
811  }
812  }
813 
814  double dsc0 = 0.5*sigUp/sigDn;
815  if(ds1 > 0.0)
816  dsc0 = 1.0 - dsc0*exp(-dsc*dsc);
817 
818  return dsc0;
819  }
820 
821 //_______________________________________________
822 double StSvtSignal::numConvInt(int nMin , int n, int numOfIntPoints, double tStep, double t)
823 {
824  double At = 0.0,lowlim = 0;
825 
826  for(int p = nMin; p <= n ; p++)
827  {
828  //lowlim = -mLowLim*(tStep/100);
829  //At = simpsonInt(mLowLim,lowlim,tStep/100,t);
830  lowlim = (p - 1)*mTimeBinSize;
831  At += simpsonInt(numOfIntPoints,lowlim,tStep/numOfIntPoints,t);
832 
833  }
834 
835  return At;
836 
837 }
838 
839 
840 //_______________________________________________
841 double StSvtSignal::analConvInt(double tim, double sigmat, double tc)
842 {
843  double a, b, c, bc5, bc4, bc3, bc2, ac1;
844  double At = 0,tb, ta, ta2,ta3, ta4,Errb,Erra,expb,expa,gaus ;
845  double sigmat2,sigmat3,sigmat4;
846 
847  sigmat2 = sigmat*sigmat;
848  sigmat3 = sigmat2*sigmat;
849  sigmat4 = sigmat3*sigmat;
850 
851  //a = 1.0/mTau_s; b = 1.0/mTau_l; c = b - a;
852  a = 90.909090909091; b = 2.0; c = -88.909090909091;
853 
854  //bc5 = b/::pow(c,5); bc4 = b/::pow(c,4); bc3 = b/(2*::pow(c,3)); bc2 = b/(6*c*c); ac1 = a/(24*c);
855 
856  bc5 = -0.00000000035999; bc4 = 0.00000003200702;
857  bc3 = -0.00000142285777; bc2 = 0.00004216833039;
858  ac1 = -0.04260395364689;
859 
860 
861  tb = tim - tc - b*sigmat2;
862  ta = tim - tc - a*sigmat2;
863  ta2 = ta*ta; ta3 = ta2*ta; ta4 = ta3*ta;
864 
865  Errb = prob2(tb,sigmat); expb = exp(-b*(tim - tc - b*(sigmat2*0.5)));
866  Erra = prob2(ta,sigmat); expa = exp(-a*(tim - tc - a*(sigmat2*0.5)));
867 
868 
869  gaus = (sigmat/::sqrt(2*M_PI))*exp(-0.5*ta2/sigmat2);
870 
871  //double c1d1 = expb*Errb;
872  //double c2d2 = expa*Erra;
873  //double c2d3 = expa*(ta*Erra + gaus);
874  //double c2d4 = expa*((ta2 + sigmat2)*Erra + ta*gaus);
875  //double c2d5 = expa*((ta3 + 3.0*ta*sigmat2)*Erra + (ta2 + 2.0*sigmat2)*gaus);
876  //double c2d6 = expa*((ta4 + 6.0*ta2*sigmat2 + 3.0*sigmat4)*Erra + (ta3 + 5.0*ta*sigmat2)*gaus);
877 
878  //At = (bc5*c1d1 - bc5*c2d2 + bc4*c2d3 - bc3*c2d4 + bc2*c2d5 - ac1*c2d6);
879 
880  At = bc5*expb*Errb + ((-bc5 + bc4*ta - bc3*(ta2 + sigmat2) + bc2*(ta3 + 3.0*ta*sigmat2)
881  - ac1*(ta4 + 6.0*ta2*sigmat2 + 3*sigmat4))*Erra + (bc4 - bc3*ta + bc2*(ta2 + 2.0*sigmat2)
882  - ac1*(ta3 + 5.0*ta*sigmat2))*gaus)*expa;
883 
884  /*
885  if(t < tc && At < 0)
886  {At = 0;
887  cout<<t<<setw(20)<<At<<endl;
888  continue;}
889  */
890 
891  //output<<t<<setw(20)<<At<<endl;
892  //cout<<t<<setw(20)<<At<<endl;
893 
894  return At;
895 }
896 //_______________________________________________
897 double StSvtSignal::simpsonInt(int numOfIntPoints,double lowlim, double step, double t)
898  {
899  double firstTerm = 0, evenTerm = 0, oddTerm = 0, lastTerm = 0, mTpr = 0;
900 
901  for(int m = 0; m <= numOfIntPoints ; m++)
902  {
903  mTpr = lowlim + m*step; //mTpr is the integration variable
904 
905  if(m == 0)
906  firstTerm = gausInput(mTpr - mTimeCenter)*pasaRes(t - mTpr); //(1/micro seconds)*(micro seconds**4)
907  else if((m != numOfIntPoints) && (m%2 == 0))
908  evenTerm += gausInput(mTpr - mTimeCenter)*pasaRes(t - mTpr);
909  else if((m != numOfIntPoints) && (m%2 == 1))
910  oddTerm += gausInput(mTpr - mTimeCenter)*pasaRes(t - mTpr);
911  else if(m == numOfIntPoints)
912  lastTerm = gausInput(mTpr - mTimeCenter)*pasaRes(t - mTpr);
913 
914  }
915 
916  double signal = (step*0.33333333333333)*(firstTerm + lastTerm + 4.0*oddTerm + 2.0*evenTerm); //(micro seconds**4)
917  //cout<<"signal = "<<signal<<endl;
918  return signal;
919 }
920 
921 //_______________________________________________
922 double StSvtSignal::gausInput(double tim)
923 {
924  double tPr2, PI,Exp1;
925 
926  PI = M_PI;
927  //cout<<"tim = "<<tim<<endl;
928  // 1.0/::sqrt(2*PI) = 0.39894228040143;
929 
930  tPr2 = 0.5*tim*tim/(mTimeWidth*mTimeWidth);
931 
932  Exp1 = (0.39894228040143/mTimeWidth)*exp(-tPr2);
933 
934  return Exp1; //[micro seconds]**(-1)
935 
936 }
937 //_______________________________________________
938 double StSvtSignal::pasaRes(double tim)
939 {
940  double a, b, c, bc5, bc4, bc3, bc2, ac1, t2, t3,t4;
941  double Fpasa = 0;
942 
943  if(tim >= 0)
944  {
945 
946  //a = 1.0/mTau_s; b = 1.0/mTau_l; c = b - a;
947  a = 90.909090909091; b = 2.0; c = -88.909090909091;
948 
949  //bc5 = b/::pow(c,5); bc4 = b/::pow(c,4); bc3 = b/(2*::pow(c,3)); bc2 = b/(6*c*c); ac1 = a/(24*c);
950  bc5 = -0.00000000035999; bc4 = 0.00000003200702;
951  bc3 = -0.00000142285777; bc2 = 0.00004216833039;
952  ac1 = -0.04260395364689;
953 
954  t2 = tim*tim;
955  t3 = tim*t2;
956  t4 = tim*t3;
957 
958  Fpasa = (bc5*exp(-b*tim)) + ((-bc5 + bc4*tim - bc3*t2 + bc2*t3 - ac1*t4)*exp(-a*tim));
959  }
960 
961  return Fpasa;
962 
963 }
964 //_______________________________________________
965 double StSvtSignal::prob1(double anOrTimeDiff , double sigma)
966 {
967  double num = 0;
968 
969  num = anOrTimeDiff/(M_SQRT2*sigma);
970 
971  double fraction = 0.5*(1 + erf(num));
972 
973  return fraction;
974 
975 }
976 //_______________________________________________
977 double StSvtSignal::freq(double num)
978 {
979  double frq = 0;
980 
981  if(num > 0.0)
982  frq = 0.5 + 0.5*erfc( num/M_SQRT2);
983  else
984  frq = 0.5* erfc(num/M_SQRT2);
985 
986  return frq;
987 
988 }
989 
990 //_______________________________________________
991 double StSvtSignal::prob2(double num , double sigma)
992 {
993  double mSum = 0, mErrf = 0;
994  double mFactorial = 0;
995  double mPowerTerm = 0;
996  double mPowerTermSquared = 0;
997  double mCountTerm = 0;
998 
999  if(fabs(num) < 5*sigma)
1000  {
1001  for(int j = 0; j<=150; j++)
1002  {
1003 
1004  if(j==0)
1005  {
1006  mFactorial = 1.0;
1007  mPowerTerm = fabs(num)/(M_SQRT2*sigma);
1008  mPowerTermSquared = mPowerTerm*mPowerTerm;
1009  mCountTerm = mPowerTerm;
1010  }
1011 
1012  else
1013  {
1014  mFactorial = (double)j*mFactorial;
1015  mPowerTerm = (mPowerTermSquared)*(double)fabs(mPowerTerm);
1016 
1017  if(j%2 == 1)
1018  mPowerTerm = (-1)*mPowerTerm;
1019 
1020  mCountTerm =(mPowerTerm /( mFactorial*(double)(2*j + 1)));
1021 
1022  }
1023 
1024  mSum += mCountTerm;
1025 
1026  }
1027 
1028  mErrf = (2.0/::sqrt(M_PI))*mSum;
1029 
1030  if(num < 0.0)
1031  mErrf = (-1.0)*mErrf;
1032  }
1033 
1034  else if(num > 5*sigma)
1035  mErrf = 1.0;
1036  else
1037  mErrf = -1.0;
1038 
1039  return 0.5*(1.0 + mErrf);
1040 }
1041 //_______________________________________________
1042 int StSvtSignal::getLowTBin()
1043 {
1044  return mLowTBin;
1045 }
1046 //_______________________________________________
1047 int StSvtSignal::getHiTBin()
1048 {
1049  return mHiTBin;
1050 }
1051 //_______________________________________________
1052 double StSvtSignal::getTimeWidth()
1053 {
1054  return mTimeWidth;
1055 }
1056 //_______________________________________________
1057 double StSvtSignal::getTimeCenter()
1058 {
1059  return mTimeCenter;
1060 }
1061 //_______________________________________________
1062 double StSvtSignal::getPeak()
1063 {
1064  return mPeakSignal*0.001; //[mV]
1065 }
1066 //_______________________________________________
1067 double StSvtSignal::getMinUnderShoot()
1068 {
1069 
1070  return mMinUnderShoot*0.001;
1071 }
1072 //_______________________________________________
1073 double StSvtSignal::getSignal(int n)
1074 {
1075  return mSignal[n]*0.001; //[mV]
1076 }
1077 //_______________________________________________
1078 void StSvtSignal::resetSignal(int lBin, int hBin)
1079 {
1080  for(int i = lBin; i <= hBin; i++)
1081  mSignal[i-1] = 0.0;
1082 }
StSvtElectronCloud
SVT electron cloud expansion routines Simulates electron cloud expansion inside of the silicon wafer...
Definition: StSvtElectronCloud.hh:18
Generated by   doxygen 1.8.5