StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StSvtElectronCloud.cc
1 /***************************************************************************
2  *
3  * $Id: StSvtElectronCloud.cc,v 1.16 2010/04/14 18:21:52 baumgart Exp $
4  *
5  * Author: Selemon Bekele
6  ***************************************************************************
7  *`
8  * Description: Calculatos cloud expansion in principal axis frame
9  *
10  ***************************************************************************
11  *
12  * $Log: StSvtElectronCloud.cc,v $
13  * Revision 1.16 2010/04/14 18:21:52 baumgart
14  * Fixed a problem for SVT hits near eta=0
15  *
16  * Revision 1.15 2009/08/10 05:21:32 baumgart
17  * Fix Minor Axis of Initial Electron Cloud Size
18  *
19  * Revision 1.14 2009/07/29 18:23:44 baumgart
20  * Modified several formulae to better account for non-infinitesimal initial hit sizes
21  *
22  * Revision 1.13 2009/06/28 04:01:46 baumgart
23  * Correction of SDD thickness
24  *
25  * Revision 1.12 2009/06/28 03:56:07 baumgart
26  * Fix of angular dependencies
27  *
28  * Revision 1.11 2009/06/11 23:17:21 baumgart
29  * Increase initial hit sizes and add projection to tSigMaj in function setInitWidths
30  *
31  * Revision 1.10 2008/11/07 15:11:31 caines
32  * Initialize variables (that REALLY do not need to be) to see if this makes valgrind happy
33  *
34  * Revision 1.9 2008/10/22 17:33:34 fine
35  * Initialize all data-members of the class StSvtElectronCloud
36  *
37  * Revision 1.8 2007/03/21 17:25:51 fisyak
38  * Ivan Kotov's drift velocities, TGeoHMatrix
39  *
40  * Revision 1.7 2005/07/23 03:37:34 perev
41  * IdTruth + Cleanup
42  *
43  * Revision 1.6 2005/02/09 14:33:35 caines
44  * New electron expansion routine
45  *
46  * Revision 1.5 2003/11/15 20:24:29 caines
47  * fixes to remove warnings at compilation
48  *
49  * Revision 1.4 2003/11/13 16:24:59 caines
50  * Further improvements to get simulator looking like reality
51  *
52  * Revision 1.1 2000/11/30 20:47:48 caines
53  * First version of Slow Simulator - S. Bekele
54  *
55  **************************************************************************/
56 #include <string.h>
57 #include "StSvtElectronCloud.hh"
58 #include "TMath.h"
59 ClassImp(StSvtElectronCloud)
60 
61 //#define EL_CL_DEBUG
62 
63 
64 StSvtElectronCloud::StSvtElectronCloud()
65  : mSigX(-1956), mSigY(-1956), mSigXY(-1956)
66  , mChargeNow(-1956), mTotCharge(-1956), mEnergy(-1956)
67  , mTheta(-1956) , mPhi(-1956), mInitPhi(-1956)
68  , mDriftVel(-1956) , mTimBinSize(-1956),mSDD_thickness(-1956) // [mm]
69  , mTrapConst(-1956), mDiffusionConst(-1956) // [mm**2/micro seconds] X and Y are calculated from this
70  , mDiffConstX(-1956) // [mm**2/micro seconds] in drift direction
71  , mDiffConstY(-1956) // [mm**2/micro seconds] in anode direction
72  , mSi_DielConst(-1956), mSi_EnergyGap(-1956)
73  , mPermitivity(-1956) // [e/(mm-V)]
74  , mSi_Mobility(-1956) // [mm**2/(V-micro seconds)]
75  , mLifeTime(-1956) // [micro seconds]
76  , mTrackId(-1956)
77 {
78  setSiliconProp(); //important - sets SVT silicon properties
79  mTrackId = 0;
80  mTotCharge = 0;
81  mChargeNow = 0;
82 
83  for(int i = 0; i < 4; i++)
84  {
85  m_dSigX[i] = 0;
86  m_dSigY[i] = 0;
87  m_dSigXY[i] = 0;
88  }
89 }
90 
91 StSvtElectronCloud::~StSvtElectronCloud()
92 {
93 
94 }
95 
96 
97 void StSvtElectronCloud::setSiliconProp()
98 {
99  mSDD_thickness = 0.28; // [mm]
100  mInitHitSize = 0.12; // [mm]
101  mLifeTime = 1000000.0; // [micro seconds]
102  mTrapConst = 0.0; // [micro seconds]
103  mDiffusionConst=0.0035; // [mm**2/micro seconds]
104  CalculateDiffXY();
105 
107  //these values in order to speed up the calculation!!!
108  mSi_DielConst = 12.0;
109  mSi_EnergyGap = 3.6; // [eV]
110  // mPermitivity = (8.854187817/1.60217733)*10000; // [e/(mm-V)]
111  mPermitivity = 55263.46959983512;
112  mSi_Mobility = 0.135; // [mm**2/(V-micro seconds)]
113 
114 }
115 
116 
117 
118 void StSvtElectronCloud::setElectronLifeTime(double tLife)
119 {
120  mLifeTime = tLife;
121 }
122 
123 
124 
125 
126 
127 void StSvtElectronCloud::CalculateDiffXY()
128 {// recalculate "diffusion" in X and Y direction after the parameters has changed
129  mDiffConstX = mDiffusionConst;
130  mDiffConstY = mDiffusionConst;
131  if(mTrapConst)
132  mDiffConstX +=mTrapConst*pow(mDriftVel,2); // [mm**2/micro seconds]
133  //cout<<"StSvtElectronCloud::CalculateDiffXY(): diffusion in anode direction = "<< mDiffConstY <<" in time direction = "<<mDiffConstX<<endl;
134 }
135 
136 void StSvtElectronCloud::setDiffusionConst(double diffConst)
137 {
138  mDiffusionConst=diffConst;
139  CalculateDiffXY();
140 }
141 
142 void StSvtElectronCloud::setDriftVelocity(double driftVel)
143 {
144  mDriftVel = driftVel; // [mm/micro seconds]
145  CalculateDiffXY();
146 }
147 
148 void StSvtElectronCloud::setTrappingConst(double trapConst)
149 {
150 
151  mTrapConst = trapConst; // [micro seconds]
152  CalculateDiffXY();
153 
154 }
155 
156 void StSvtElectronCloud::setPar(double energy,double theta, double phi, double timBinSize,int trackId)
157 {
158  mTrackId = trackId;
159  mTimBinSize = timBinSize;
160  mEnergy = energy;
161  mTheta = theta;
162  mInitPhi = phi;
163  //mTheta = 2*acos(-1.0)/360;
164 
165  //mTotCharge = mEnergy/mSi_EnergyGap;
166  mTotCharge = mEnergy*0.27777777777777; // in number of electrons, ~25000 electrons for MIPs
167  mTotCharge = mTotCharge*(1- 2*mInitHitSize*fabs(sin(mTheta))/(3*mSDD_thickness)); // Fix eta-dependence
168  //cout<<"mTotCharge = "<<mTotCharge<<endl;
169  }
170 
171 void StSvtElectronCloud::setInitWidths(double w1, double w2)
172 {
173  // cout<<"======= NEW HIT========="<<endl;
174  double tSigMaj,tSigMin;
175  mSigX=0; mSigY=0; mSigXY=0;
176 
177  //tSigMaj = 0.288675134*fabs(mSDD_thickness*tan(mTheta)); // [mm]
178  tSigMaj = 0.288675134*(fabs(mSDD_thickness*tan(mTheta))+2*w1*cos(mTheta)); // [mm]
179 
180  if (tSigMaj<w1) { //almost perpendicular
181  tSigMaj = w1; //initial size cannot be smaller than minimal
182  mTheta=0.;
183  mPhi = 0.;
184  }
185  tSigMin= w1;
186 
187  if ((1. - tSigMin/tSigMaj) < 0.001 || (fabs(mPhi)< 0.001)) mPhi=0.;
188 
189  //calculate initial size in XY
190  if (mPhi==0.){
191  mSigX = tSigMaj;
192  mSigY = tSigMin;
193  mSigXY = 0.;
194  }
195  else{
196  double C2=cos(mPhi);
197  double S2=sin(mPhi);
198  mSigX = sqrt(fabs(tSigMaj*tSigMaj*C2*C2 + tSigMin*tSigMin*S2*S2)); //fabs --numerical safety
199  mSigY = sqrt(fabs(tSigMaj*tSigMaj*S2*S2 + tSigMin*tSigMin*C2*C2));
200  mSigXY =sqrt(fabs(C2*S2*(tSigMaj*tSigMaj-tSigMin*tSigMin)));
201  if(((1. - tSigMin/tSigMaj) < 0.001) || (fabs(mPhi)< 0.001)){
202  mPhi=0.;
203  mSigXY = 0.;
204  }
205  }
206  /*
207  cout<<"major="<<tSigMaj<<endl;
208  cout<<"initial mPhi = "<<mPhi<<endl;
209  cout<<"initial width along drift= "<<sqrt(sigmaXSq)<<endl;
210  cout<<"initial width along anode = "<<sqrt(sigmaYSq)<<endl;
211  cout<<"initial XY = "<<sqrt(sigmaXYSq)<<endl;
212  */
213 
214  for(int i = 0; i < 4; i++)
215  {
216  m_dSigX[i] = 0;
217  m_dSigY[i] = 0;
218  m_dSigXY[i] = 0;
219  }
220 
221  GetDerivatives( m_dSigX[0],m_dSigY[0], m_dSigXY[0],mSigX,mSigY,mSigXY,0.);
222  }
223 
224 void StSvtElectronCloud::CalcExpansion(double mTc)
225 {//mTc is in timebins
226  //cout<<"#### NEW HIT####"<<endl;
227  //cout <<"mDiffConstX="<<mDiffConstX <<", mDiffConstY="<<mDiffConstY <<endl;
228  const int AdBashDiv=30; //number od steps ber timebin for adams-bashfort
229  const int RungeDiv =200; //number of steps inside one AdBashStep
230 
231  //the phi has to be from -2pi to 2pi
232  //now remember in which quandrant was phi
233  double sign = (mInitPhi<0.)?-1.:1.;
234  mPhi=fabs(mInitPhi);
235  int quad=(int)floor(mPhi/TMath::PiOver2());
236  //cout<<"quandrant="<<quad<<endl;
237  switch (quad){
238  case 0:break;
239  case 1:
240  mPhi=TMath::Pi()-mPhi;
241  break;
242  case 2:
243  mPhi=mPhi-TMath::Pi();
244  break;
245  case 3:
246  mPhi=TMath::TwoPi()-mPhi;
247  break;
248  default: //just in case
249  mPhi=asin(sin(mPhi));
250  break;
251  }
252 
253 
254  //set initial sizes
255  mPhi = mInitPhi;
256 
257  //setInitWidths(0.002, 0.002);
258  setInitWidths(mInitHitSize,mInitHitSize); // Stephen tune to Run 7 Au+Au
259 
260  double steps=mTc*AdBashDiv*RungeDiv;
261  int isteps=(int)floor(steps + 0.5);
262  /*
263  cout<<" time="<<mTc<<" isteps="<<isteps<<endl;
264  cout<<" ..initial size: mSigX="<<mSigX<<" mSigY="<<mSigY<<" mSigXY="<<mSigXY<<endl;
265  cout<<" ..comparing : SigT="<<getSigmaDrift()<<" Anode="<<getSigmaAnode()<<endl;
266  cout<<" ..comparing : Major="<<getSigmaMajor()<<" Minor="<<getSigmaMinor()<<" Phi="<<getPhi()<<endl;
267  */
268  if (isteps<=3*RungeDiv){//very short drift time -handle all by runge-kutta
269  //cout<<" going for short runge-kutta: steplen="<<(mTimBinSize/((double)AdBashDiv*(double)RungeDiv))<<endl;
270  runge_kutta4(isteps,0., mTimBinSize/((double)AdBashDiv*(double)RungeDiv),0);
271  return;
272  }
273 
274  //now four runge-kutta runs to get input for adams-bashfort
275  //first runs to whole number of adams-basfort steps
276  int tAdBashSteps=isteps/RungeDiv;
277  int firstNumStep=isteps-tAdBashSteps*RungeDiv;
278  tAdBashSteps-=3;
279  double steplen=mTimBinSize/((double)AdBashDiv*(double)RungeDiv);
280 
281 #ifdef EL_CL_DEBUG
282  cout<<" runge-kutta: firstNumSteps="<<firstNumStep<<" steplen="<<steplen<<endl;
283 #endif
284  runge_kutta4(firstNumStep, 0.,steplen,0);
285  runge_kutta4(RungeDiv,steplen*(double)firstNumStep ,steplen,1);
286  runge_kutta4(RungeDiv, steplen*((double)firstNumStep+(double)RungeDiv), steplen,2);
287  runge_kutta4(RungeDiv, steplen*((double)firstNumStep+2.*(double)RungeDiv), steplen,3);
288 
289  //Now run Adams-Bashfort for the rest
290 #ifdef EL_CL_DEBUG
291  cout<<" going for adamsBushFort: steps:="<<tAdBashSteps<<endl;
292 #endif
293  double Adsteplen=mTimBinSize/(double)AdBashDiv;
294 #ifdef EL_CL_DEBUG
295  cout<<"T0="<<steplen*((double)firstNumStep+3.*(double)RungeDiv)<<endl;
296 #endif
297  adamsBushFort(tAdBashSteps,steplen*((double)firstNumStep+3.*(double)RungeDiv),Adsteplen);
298 
299  //return to the propper quadrant
300  switch (quad){
301  case 0:break;
302  case 1:
303  mPhi=TMath::Pi()-mPhi;
304  break;
305  case 2:
306  mPhi=mPhi+TMath::Pi();
307  break;
308  case 3:
309  mPhi=TMath::TwoPi()-mPhi;
310  break;
311  default:
312  break;
313  }
314 
315  mPhi=sign*mPhi;
316 }
317 
318 void StSvtElectronCloud::GetDerivatives(double &dSx,double &dSy,double &dSxy,double SX,double SY,double SXY,double time)
319 { //these are derivatives for sigmas^3 - it's needed for numerical stability
320 
321  // double sr, sr2, sr3,sr4, sd, sd2, cons,cons3;
322  //double fun1 ,multFactor, denominator,numerator;
323 #ifdef EL_CL_DEBUG
324  cout<<"SX, SY, SXY : "<<SX<<", "<<SY<<", "<<SXY<<endl;
325 #endif
326 
327  //----first get major, minor and angle
328  double R=SX*SX + SY*SY;
329  double S=sqrt(pow((SX*SX - SY*SY),2) + 4.*pow(SXY*SXY,2));
330  double sigMaj = sqrt(0.5*fabs(R + S));
331  double sigMin = sqrt(0.5*fabs(R - S)); //safety fabs
332  double phi;
333  if ((1. - sigMin/sigMaj) < 0.001)phi=0;
334  else
335  {
336  double ar=(SX*SX-SY*SY)/(sigMaj*sigMaj-sigMin*sigMin);
337  #ifdef EL_CL_DEBUG
338  cout<<" ar="<<(double)ar<<endl;
339 #endif
340  if (ar>.999999) ar=(double)1e0;
341  if (ar<-.999999) ar=(double)-1e0;
342  phi=0.5*acos(ar);
343  }
344 #ifdef EL_CL_DEBUG
345  cout<<" R S sigMaj sigMin phi: "<<S<<", "<< sigMaj<<", "<< sigMin<<", "<< phi<<endl;
346 #endif
347  //----calc derivatives
348  const double K1 = 0.43;
349  const double K2 = 0.25;
350  const double K3 = 0.58;
351 
352  mChargeNow = mTotCharge*exp(-time/mLifeTime);
353  //multFactor = (1.0/(4*pi))*(mSi_Mobility/(mSi_DielConst*mPermitivity));
354  double multFactor = 0.00000001619961;
355 
356 
357  double vd=sigMin/mSDD_thickness;
358  double r = sigMin/sigMaj;
359  double r2=r*r;
360 
361  double tmp=K1*vd/K2;
362  tmp=1. + pow(tmp,1.5);
363  double denominator=pow(tmp,2./3.);
364 
365  tmp=(1.+1./(vd*vd));
366  double numerator = K1*pow(r,2./3.) - (K3/4.)*log(r2*r2 + 1./(tmp*tmp));
367 #ifdef EL_CL_DEBUG
368  cout<<" vd="<<vd<<" r2="<<r2<<" den="<<denominator<<" num1="<<numerator;
369 #endif
370  double f1 = multFactor*mChargeNow*(numerator/denominator);
371 
372  //now func2
373  numerator = K3 - (K3-K1)*sqrt(r);
374  double f2 = multFactor*mChargeNow*(numerator/denominator);
375 #ifdef EL_CL_DEBUG
376  cout<<" num2="<<numerator<<endl;
377 #endif
378  double sumFunc=f1+f2;
379  double diffFunc=f1-f2;
380  double C=cos(2.*phi);
381 
382  dSx = 1.5*SX*(2*mDiffConstX + 0.5*(sumFunc + diffFunc*C));
383  dSy = 1.5*SY*(2*mDiffConstY + 0.5*(sumFunc - diffFunc*C));
384  dSxy = 1.5*SXY*(2.*(sin(2*phi)/C)*(mDiffConstX*pow(sin(phi),2) - mDiffConstY*pow(cos(phi),2))
385  + 0.5*diffFunc*sin(2.*phi));
386 #ifdef EL_CL_DEBUG
387  cout<<" dSx,dSy,dSxy= "<<dSx<<","<<dSy<<","<<dSxy<<","<<endl;
388 #endif
389 
390  /*-exponential test
391  dSx =3*SX*SX*SX;
392  dSy = 3*SY*SY*SY;
393  dSxy=3*SXY*SXY*SXY;
394  */
395 }
396 
397 void StSvtElectronCloud::runge_kutta4(int steps, double t0, double steplen,int save)
398 {
399  //runge kutta 4th order
400  double m1 = 0.0, m2 = 0.0, m3 = 0.0, m4 = 0.0;
401  double n1 = 0.0, n2 = 0.0, n3 = 0.0, n4 = 0.0;
402  double o1 = 0.0, o2 = 0.0, o3 = 0.0, o4 = 0.0;
403 
404  double tim = t0;
405 
406 
407 
408  /*
409  cout<<"####Runge-kutta -entry: tim = "<<tim<<" #############"<<endl;
410  cout<<"step length = "<<steplen<<endl;
411  cout<<"theta = "<<mTheta<<" phi = "<<phi<<endl;
412  cout<<"sigma1Sq = "<<sigma1Sq<<" sigma2Sq = "<<sigma2Sq<<endl;
413  cout<<" width along drift= "<<sqrt(mSigX)<<" width along anode = "<<sqrt(mSigY)<<endl;
414  cout<<"initial XY = "<<sqrt(mSigXY)<<endl;
415  */
416 
417  //phi for first loop
418  const double third=1./3.;
419  if (steps>0)
420  for(int i = 1; i <= steps; i++)
421  {
422  /*
423  cout<<" +++++ Runge Kutta -step"<<i ;
424  cout<<" SigT="<<getSigmaDrift()<<" Anode="<<getSigmaAnode()<<endl;
425  cout<<" Major="<<getSigmaMajor()<<" Minor="<<getSigmaMinor()<<" Phi="<<getPhi()<<endl;
426  */
427  double xx=mSigX*mSigX*mSigX;
428  double yy=mSigY*mSigY*mSigY;
429  double xy=mSigXY*mSigXY*mSigXY;
430 
431  GetDerivatives(m1,n1,o1,mSigX,mSigY,mSigXY,tim);
432  //cout<<" dXX="<< m1<<" dYY="<<n1 <<" dXY="<< o1<<endl;
433 
434  double tmpX=pow(fabs(xx+ 0.5*steplen*m1),third);
435  double tmpY=pow(fabs(yy+ 0.5*steplen*n1),third);
436  double tmpXY=pow(fabs(xy+ 0.5*steplen*o1),third);
437  GetDerivatives(m2,n2,o2,tmpX,tmpY,tmpXY,tim+0.5*steplen);
438 
439  tmpX=pow(fabs(xx+0.5*steplen*m2),third);
440  tmpY=pow(fabs(yy+ 0.5*steplen*n2),third);
441  tmpXY=pow(fabs(xy+ 0.5*steplen*o2),third);
442  GetDerivatives(m3,n3,o3,tmpX,tmpY,tmpXY,tim+0.5*steplen);
443 
444  tmpX=pow(fabs(xx+steplen*m3),third);
445  tmpY=pow(fabs(yy+ steplen*n3),third);
446  tmpXY=pow(fabs(xy+ steplen*o3),third);
447  GetDerivatives(m4,n4,o4,tmpX,tmpY,tmpXY,tim+steplen);
448 
449 
450  mSigX = pow(fabs(xx + ((steplen/6.)*(m1 + 2.*(m2 + m3) + m4))),third);
451  mSigY = pow(fabs(yy + ((steplen/6.)*(n1 + 2.*(n2 + n3) + n4))),third);
452  mSigXY = pow(fabs(xy + ((steplen/6.)*(o1 + 2.*(o2 + o3) + o4))),third);
453 
454 
455 
456 
457  tim = tim + steplen;
458  }//for loop
459 
460  if ((mPhi==0.)||(mPhi==TMath::PiOver2()))mSigXY = 0.;
461  else
462  {
463  double R=mSigX*mSigX + mSigY*mSigY;
464  double S=sqrt(pow((mSigX*mSigX - mSigY*mSigY),2) + 4.*pow(mSigXY*mSigXY,2));
465  double sigMaj2 = 0.5*fabs(R + S);
466  double sigMin2 = 0.5*fabs(R - S); //safety fabs
467  if (((1. - sigMin2/sigMaj2) < 0.00001)|| (fabs(mPhi)< 0.001))mPhi=0;
468  else{
469  double ar=(mSigX*mSigX-mSigY*mSigY)/(sigMaj2-sigMin2);
470  if (ar>.999999) ar=(double)1e0; //due to numerical problems
471  if (ar<-.999999) ar=(double)-1e0;
472  mPhi=0.5*acos(ar);
473  }
474  }
475 
476  /*
477  cout<<"---leaving Runge-kutta----"<<endl;
478  cout<<"sigma drift = "<<sqrt(mSigX)<<endl;
479  cout<<"sigma anode= "<<sqrt(mSigY)<<endl;
480  cout<<"Sigma XY = "<<sqrt(mSigXY)<<endl;
481  cout<<"phi ="<<mPhi<<endl;
482  cout<<"###########################"<<endl;
483  */
484 
485  GetDerivatives(m1,n1,o1,mSigX,mSigY,mSigXY,tim);
486  m_dSigX[save] = m1;
487  m_dSigY[save] = n1;
488  m_dSigXY[save] = o1;
489 
490  mChargeNow = mTotCharge*exp(-tim/mLifeTime);
491 }
492 
493 void StSvtElectronCloud::adamsBushFort(int steps, double t0, double steplen)
494 {
495  double tim = t0;
496  const double third=1./3.;
497 
498  while (steps>0)
499  {
500 #ifdef EL_CL_DEBUG
501  cout<<"time ="<<tim<<endl;
502 #endif
503  double xx=mSigX*mSigX*mSigX;
504  double yy=mSigY*mSigY*mSigY;
505  double xy=mSigXY*mSigXY*mSigXY;
506 #ifdef EL_CL_DEBUG
507  cout<<"sigX="<<mSigX<<" sigY="<<mSigY<<" sigXY="<<mSigXY<<endl;
508 #endif
509 
510  double preX = pow(fabs(xx + (steplen/24.)*(-9.0*m_dSigX[0] + 37.0*m_dSigX[1] - 59.0*m_dSigX[2] + 55.0*m_dSigX[3])),third);
511  double preY = pow(fabs(yy + (steplen/24.)*(-9.0*m_dSigY[0] + 37.0*m_dSigY[1] - 59.0*m_dSigY[2] + 55.0*m_dSigY[3])),third);
512  double preXY = pow(fabs(xy + (steplen/24.)*(-9.0*m_dSigXY[0] + 37.0*m_dSigXY[1] - 59.0*m_dSigXY[2] + 55.0*m_dSigXY[3])),third);
513 #ifdef EL_CL_DEBUG
514  cout<<"preX="<<preX<<" preY="<<preY<<" preXY="<<preXY<<endl;
515 #endif
516 
517  tim = tim + steplen;
518  double dX=0;
519  double dY=0;
520  double dXY=0;
521  GetDerivatives(dX,dY,dXY,preX,preY,preXY,tim);
522 
523  //get the corrector...also final value
524  mSigX=pow(fabs( xx + (steplen/24.)*(m_dSigX[1]-5.0*m_dSigX[2] + 19.0*m_dSigX[3] + 9.0*dX)),third);
525  mSigY=pow(fabs( yy + (steplen/24.)*(m_dSigY[1]-5.0*m_dSigY[2] + 19.0*m_dSigY[3] + 9.0*dY)),third);
526  mSigXY=pow(fabs( xy + (steplen/24.)*(m_dSigXY[1]-5.0*m_dSigXY[2] + 19.0*m_dSigXY[3] + 9.0*dXY)),third);
527 
528 
529 
530  GetDerivatives(dX,dY,dXY,mSigX,mSigY,mSigXY,tim);
531  //next CE
532 
533  mSigX=pow(fabs( xx + (steplen/24.)*(m_dSigX[1]-5.0*m_dSigX[2] + 19.0*m_dSigX[3] + 9.0*dX)),third);
534  mSigY=pow(fabs( yy + (steplen/24.)*(m_dSigY[1]-5.0*m_dSigY[2] + 19.0*m_dSigY[3] + 9.0*dY)),third);
535  mSigXY=pow(fabs( xy + (steplen/24.)*(m_dSigXY[1]-5.0*m_dSigXY[2] + 19.0*m_dSigXY[3] + 9.0*dXY)),third);
536 
537 
538 
539  GetDerivatives(dX,dY,dXY,mSigX,mSigY,mSigXY,tim);
540 
541  //shift the values
542  m_dSigX[3] = dX;
543  m_dSigY[3] = dY;
544  m_dSigXY[3] = dXY;
545  for(int i = 0; i < 3; i++){
546  m_dSigX[i] = m_dSigX[i+1];
547  m_dSigY[i] = m_dSigY[i+1];
548  m_dSigXY[i] = m_dSigXY[i+1];
549  }
550 
551  mChargeNow = mTotCharge*exp(-tim/mLifeTime);
552  steps--;
553  }
554 
555  //get phi
556  if (mPhi==0.||(mPhi==TMath::PiOver2())) mSigXY = 0.;
557  else
558  {
559  double R=mSigX*mSigX + mSigY*mSigY;
560  double S=sqrt(pow((mSigX*mSigX - mSigY*mSigY),2) + 4.*pow(mSigXY*mSigXY,2));
561  double sigMaj2 = 0.5*fabs(R + S);
562  double sigMin2 = 0.5*fabs(R - S); //safety fabs
563  if (((1. - sigMin2/sigMaj2) < 0.00001)|| (fabs(mPhi)< 0.001))mPhi=0;
564  else{
565  double ar=(mSigX*mSigX-mSigY*mSigY)/(sigMaj2-sigMin2);
566  if (ar>.999999) ar=(double)1e0;
567  if (ar<-.999999) ar=(double)-1e0;
568  mPhi=0.5*acos(ar);
569  }
570  }
571 }
572 
573 
574 
575 double StSvtElectronCloud::getPhi()
576 {
577  return mPhi;
578 }
579 
580 double StSvtElectronCloud::getChargeAtAnode()
581 {
582  //cout<<" mChargeNow = "<<mChargeNow<<endl;
583 return mChargeNow;
584 
585 }
586 
587 
588 double StSvtElectronCloud::getSigmaDrift()
589 {
590  return mSigX;
591 }
592 
593 double StSvtElectronCloud::getSigmaAnode()
594 {
595  return mSigY;
596 }
597 
598 double StSvtElectronCloud::getSigmaCorr()
599 {
600  return mSigXY;
601 }
602 
603 double StSvtElectronCloud::getSigmaMajor()
604 {
605  double c2f=mSigX*mSigX + mSigY*mSigY;
606  double s2f=sqrt(pow((mSigX*mSigX - mSigY*mSigY),2) + 4.*pow(mSigXY*mSigXY,2));
607  return sqrt(0.5*fabs(c2f+s2f));
608 }
609 
610 double StSvtElectronCloud::getSigmaMinor()
611 {
612  //return sqrt(mSigY);
613  return sqrt(0.5*fabs(((mSigX*mSigX + mSigY*mSigY) - sqrt(pow((mSigX*mSigX - mSigY*mSigY),2) + 4.*pow(mSigXY*mSigXY,2)))));
614 }
615 
616 
617 #ifdef EL_CL_DEBUG
618 #undef EL_CL_DEBUG
619 #endif
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