StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFtpcLaserTrafo.cc
1 // $Id: StFtpcLaserTrafo.cc,v 1.8 2010/01/14 18:44:28 jcs Exp $
2 //
3 // $Log: StFtpcLaserTrafo.cc,v $
4 // Revision 1.8 2010/01/14 18:44:28 jcs
5 // change padtrans so that the code is the same as StFtpcClusterFinder::padtrans
6 //
7 // Revision 1.7 2009/12/09 14:40:01 jcs
8 // update comment
9 //
10 // Revision 1.6 2007/01/22 13:08:15 jcs
11 // replace cout with LOG
12 //
13 // Revision 1.5 2007/01/19 08:23:00 jcs
14 // replace true, false with kTRUE, kFALSE
15 //
16 // Revision 1.4 2007/01/16 15:04:17 jcs
17 // reolace "return FALSE" with "return false" as temporary fix so that dev compiles
18 //
19 // Revision 1.3 2006/04/12 08:46:01 jcs
20 // initialize p0,p1,p2
21 //
22 // Revision 1.2 2006/03/15 15:13:57 jcs
23 // add lines for listing CVS update info
24 //
25 
26 #include "StFtpcLaserTrafo.hh"
27 #include "TSystem.h"
28 #include "StMessMgr.h"
29 #include "TMath.h"
30 #include "PhysicalConstants.h"
31 
32 const float normgas = 50.0;
33 const Float_t Grad2Rad=1/(180/TMath::Pi());
34 //---------------------------------------------------------------
35 
36 StFtpcLaserTrafo::StFtpcLaserTrafo()
37 {
38  // default constructor !!!
39 }
40 //---------------------------------------------------------------
41 
42 StFtpcLaserTrafo::StFtpcLaserTrafo(StFtpcDbReader *getdb,StFtpcParamReader *getparam,float getdeltat0,float getdeltagas, float getmicropertime, float getdeltap, float getmbfield, float getTZero)
43 {
44  LOG_INFO<<"StFtpcLaserTrafo..."<<endm;
45  mDb=getdb; mParam=getparam;
46 
47  mBField = getmbfield;
48 
49  deltat0=getdeltat0;
50  deltagas=getdeltagas;
51  tZero = getTZero;
52 
53  LOG_INFO<<"d_t0 = "<<deltat0<<" | d_gas(Ar) = "<<(float) deltagas<<endm;
54  micropertime = getmicropertime;
55  LOG_INFO<<"Microseconds per timebin = "<<micropertime<<endm;
56  LOG_INFO<<"==============================="<<endm;
57  LOG_INFO<<" "<<endm;
58 
59  deltap = getdeltap;
60 
61  pradius = new Double_t[mParam->numberOfDriftSteps()
62  *mDb->numberOfPadrowsPerSide()];
63  pdeflection = new Double_t[mParam->numberOfDriftSteps()
64  *mDb->numberOfPadrowsPerSide()];
65 }
66 
67 //---------------------------------------------------------------
68 
69 StFtpcLaserTrafo::~StFtpcLaserTrafo()
70 {
71  delete[] pradius; // release the pradius array
72  delete[] pdeflection; // release the pdeflection array
73  //LOG_DEBUG<<"StFtpcLaserTrafo deconstructed"<<endm;
74 }
75 
76 //---------------------------------------------------------------
77 
78 //double StFtpcLaserTrafo::vd_gas(float rad, float gasmix)
79 //{
80  //return vd_gas_slope(rad)*(gasmix-normgas)+vd_gas_y(rad);
81  //return vd_gas_slope(rad)*(gasmix-normgas);
82 //}
83 
84 //---------------------------------------------------------------
85 
86 double StFtpcLaserTrafo::vd_gas(float rad)
87 {
88  //return vd_gas_slope(rad)*(gasmix-normgas)+vd_gas_y(rad);
89  return vd_gas_slope(rad)*deltagas;
90 }
91 
92 //---------------------------------------------------------------
93 
94 double StFtpcLaserTrafo::vd_gas_slope(float rad)
95 {
96  double p0 = 0.0;
97  double p1 = 0.0;
98  double p2 = 0.0;
99 
100  if ( mBField == 0) {
101  p0=0.8327;
102  p1=1.613;
103  p2=0.00378;
104  }
105  if ( abs(mBField) == 1) {
106  p0=0.7775;
107  p1=1.598;
108  p2=0.003616;
109  }
110 
111  return (p0*1/TMath::Power(rad,p1)+p2);
112 }
113 
114 //---------------------------------------------------------------
115 // B=1
116 double StFtpcLaserTrafo::vd_gas_y(float rad)
117 {
118  double p0=13.05;
119  double p1=1.08;
120 
121  return p0*1/TMath::Power(rad,p1);
122 }
123 
124 //---------------------------------------------------------------
125 
126 //double StFtpcLaserTrafo::lor_gas(float rad, float gasmix)
127 //{
128  //return lor_gas_slope(rad)*(gasmix-normgas)+lor_gas_y(rad);
129  //return lor_gas_slope(rad)*(gasmix-normgas);
130 //}
131 
132 //---------------------------------------------------------------
133 
134 double StFtpcLaserTrafo::lor_gas(float rad)
135 {
136  //return lor_gas_slope(rad)*(gasmix-normgas)+lor_gas_y(rad);
137  return lor_gas_slope(rad)*deltagas;
138 }
139 
140 //---------------------------------------------------------------
141 
142 double StFtpcLaserTrafo::lor_gas_slope(float rad)
143 {
144  double p0 = 0.0;
145  double p1 = 0.0;
146  double p2 = 0.0;
147 
148  if ( mBField == 0 ) {
149  p0=0.0001995;
150  p1=0.1211;
151  p2=-0.0001324;
152  }
153  if ( abs(mBField) == 1) {
154  p0=8.916;
155  p1=2.723;
156  p2=0.0775;
157  }
158 
159  return (p0*1/TMath::Power(rad,p1)+p2);
160 }
161 
162 
163 //---------------------------------------------------------------
164 // B=1
165 double StFtpcLaserTrafo::lor_gas_y(float rad)
166 {
167 
168  // nicht besonders gut der fit !!!)
169 
170  double p0=8.005;
171  double p1=-0.9733;
172  double p2=0.09855;
173  double p3=-0.005031;
174  double p4=0.0001277;
175  double p5=-1.281/1000000;
176 
177  return p0+p1*rad+p2*TMath::Power(rad,2)+p3*TMath::Power(rad,3)+p4*TMath::Power(rad,4)+p5*TMath::Power(rad,5);
178 }
179 //---------------------------------------------------------------
180 
181 // include Gas variation !!!
182 
183 int StFtpcLaserTrafo::calcpadtrans()
184 {
185  int i, j, v_buf, padrow;
186  double t_last, t_next, r_last, r_next, e_now, v_now, psi_now;
187  double step_size;
188 
189  LOG_INFO<<"calcpadtrans ..."<<endm;
190 
191  step_size=((float) mDb->numberOfTimebins()
192  / (float) mParam->numberOfDriftSteps());
193 
194 
195  LOG_INFO<<"Pressure & temperature correction = "<<deltap<<endm;
196  LOG_INFO<<"======================================================"<<endm;
197 
198 
199  LOG_INFO<<"RadiusTimesField = "<<mDb->radiusTimesField()<<endm;
200  LOG_INFO<<"vdrift(0,0) = "<<mDb->magboltzVDrift(0,0)<<endm;
201  LOG_INFO<<"lorangle(0,0) = "<<mDb->magboltzDeflection(0,0)<<endm;
202  LOG_INFO<<" "<<endm;
203 
204  for (padrow=0; padrow<mDb->numberOfPadrowsPerSide(); padrow++)
205  {
206  /* determine starting values */
207  t_last=0;
208  v_buf=0;
209  r_last=mDb->sensitiveVolumeOuterRadius();
210  pradius[padrow]=mDb->sensitiveVolumeOuterRadius();
211  pdeflection[padrow]=0;
212  e_now = mDb->radiusTimesField() / (0.5*r_last);
213  for(j=v_buf; mDb->magboltzEField(j) < e_now
214  && j<(mDb->numberOfMagboltzBins()-1); j++);
215  if(j<1 || j>mDb->numberOfMagboltzBins())
216  {
217  gMessMgr->Message("", "E", "OST") << "Error 1: j=" << j << ", v_buf=" << v_buf << " e_drift=" << mDb->magboltzEField(j) << ", e_now=" << e_now << endm;
218  return kFALSE;
219  }
220  v_buf=j-1;
221  v_now=((mDb->magboltzVDrift(v_buf, padrow)
222  +deltap*mDb->magboltzdVDriftdP(v_buf, padrow))
223  *(mDb->magboltzEField(j)-e_now)
224  +(mDb->magboltzVDrift(j, padrow)
225  +deltap*mDb->magboltzdVDriftdP(j, padrow))
226  *(e_now-mDb->magboltzEField(v_buf)))
227  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
228  psi_now=((mDb->magboltzDeflection(v_buf,padrow)
229  +deltap*mDb->magboltzdDeflectiondP(v_buf,padrow))
230  *mParam->lorentzAngleFactor()
231  *(mDb->magboltzEField(j)-e_now)
232  +(mDb->magboltzDeflection(j,padrow)
233  +deltap*mDb->magboltzdDeflectiondP(j,padrow))
234  *mParam->lorentzAngleFactor()
235  *(e_now-mDb->magboltzEField(v_buf)))
236  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
237 
238  //LOG_DEBUG<<psi_now<<endm;
239 
240  // gas correction (only use by fullfield !)
241  if (deltagas!=0)
242  {
243  //LOG_DEBUG<<"v_now before "<<v_now<<endm;
244  v_now=v_now+vd_gas(r_last);
245  //LOG_DEBUG<<"v_now after "<<v_now<<endm;
246  //LOG_DEBUG<<"pnow "<<psi_now<<endm;
247  //LOG_DEBUG<<Grad2Rad*lor_gas(r_last)<<" "<<lor_gas(r_last)<<endm;
248  //psi_now=psi_now+lor_gas(r_last);
249  // Is 1st approximation reversed FullField correct !????
250  //psi_now=-(psi_now+lor_gas(r_last));
251 
252  if (mBField == -1 )
253  psi_now=psi_now-lor_gas(r_last); // ???? B=-1 !!!
254  else if (mBField == 0 || mBField == 1 )
255  psi_now=psi_now+lor_gas(r_last); // B=0,B=+1
256  else {
257  LOG_ERROR<<"psi_now not defined for mBField = "<<mBField<<endm;
258  }
259 
260  //LOG_DEBUG<<"psi_now "<<psi_now<<endm;
261  //v_now=v_now*cos(Grad2Rad*psi_now);
262 
263  //***************************
264  //* ??? so 1. Naeherung ??? *
265  //***************************
266  //LOG_DEBUG<<"v_now "<<v_now<<endm;
267  }
268 
269  for (i=0; i<mParam->numberOfDriftSteps()
270  && e_now < mDb->magboltzEField(mDb->numberOfMagboltzBins()-2)
271  ; i++)
272  {
273  t_next = t_last + step_size;
274  /* first guess for r_next: */
275  //r_next = r_last - v_now * step_size * mDb->microsecondsPerTimebin();
276  r_next = r_last - v_now * step_size * micropertime;
277  e_now = mDb->radiusTimesField() / (0.5*(r_last+r_next));
278 
279  for(j=v_buf; mDb->magboltzEField(j) < e_now
280  && (j<mDb->numberOfMagboltzBins()-1); j++);
281 
282  if(j<1 || j>mDb->numberOfMagboltzBins())
283  {
284  gMessMgr->Message("", "E", "OST") << "Error 2: j=" << j << ", v_buf=" << v_buf << " e_drift=" << mDb->magboltzEField(j) << ", e_now=" << e_now << endm;
285  return kFALSE;
286  }
287 
288  v_buf=j-1;
289  v_now=((mDb->magboltzVDrift(v_buf, padrow)
290  +deltap*mDb->magboltzdVDriftdP(v_buf, padrow))
291  *(mDb->magboltzEField(j)-e_now)
292  +(mDb->magboltzVDrift(j, padrow)
293  +deltap*mDb->magboltzdVDriftdP(j, padrow))
294  *(e_now-mDb->magboltzEField(v_buf)))
295  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
296  psi_now=((mDb->magboltzDeflection(v_buf,padrow)
297  +deltap*mDb->magboltzdDeflectiondP(v_buf,padrow))
298  *mParam->lorentzAngleFactor()
299  *(mDb->magboltzEField(j)-e_now)
300  +(mDb->magboltzDeflection(j,padrow)
301  +deltap*mDb->magboltzdDeflectiondP(j,padrow))
302  *mParam->lorentzAngleFactor()
303  *(e_now-mDb->magboltzEField(v_buf)))
304  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
305 
306  // gas correction (only use by fullfield !)
307  if (deltagas!=0)
308  {
309 
310  v_now=v_now+vd_gas(r_last); // aenderung von rlast -> rnext !!!???
311  //psi_now=psi_now+lor_gas(r_last);
312  // 1. Neherung reversed FullField richtig !????
313  //psi_now=-(psi_now+lor_gas(r_last));
314 
315  //psi_now=psi_now-lor_gas(r_last); // ???? B=-1 !!!
316 
317  psi_now=psi_now+lor_gas(r_last); // B=0,B=+1
318  //v_now=v_now*cos(Grad2Rad*psi_now);
319  //LOG_DEBUG<<(lor_gas(r_next)-lor_gar_nexts(r_last))*180/TMath::Pi()<<endm;
320  //LOG_DEBUG<<"r_last = "<<r_last<<" | r_next = "<<r_next<<endm;
321  }
322 
323  /* correct r_next: */
324  //r_next = r_last - v_now * step_size *mDb->microsecondsPerTimebin();
325  //LOG_DEBUG<<"r_last = "<<r_last<<" | r_next = "<<r_next<<" | r_next_new = "<< r_last - v_now * step_size * micropertime<<" | psi = "<<psi_now<<endm;
326  r_next = r_last - v_now * step_size * micropertime;
327  pradius[padrow+mDb->numberOfPadrowsPerSide()*(i+1)]=r_next;
328  //LOG_DEBUG<<padrow+mDb->numberOfPadrowsPerSide()*(i+1)<<endl;
329  pdeflection[padrow+mDb->numberOfPadrowsPerSide()*(i+1)]
330  =pdeflection[padrow+mDb->numberOfPadrowsPerSide()*i]
331  +((r_last-r_next)*tan(degree * psi_now)/r_last);
332 
333  t_last=t_next;
334  r_last=r_next;
335  }
336  }
337  return kTRUE;
338 }
339 
340 //---------------------------------------------------------------
341 
342 int StFtpcLaserTrafo::padtrans(int iRow,int iSec,float timepos, float padpos,float *x1,float *y1)
343 {
344  int PadtransPerTimebin;
345  int PadtransLower;
346  float PhiDeflect, TimeCoordinate;
347  iRow = iRow - 1;
348  iSec = iSec - 1;
349 
350  TimeCoordinate = timepos + 0.5; /*time start at beginning of bin 0*/;
351 
352  // Laser t0 from Calibrations_ftpc/ftpcElectronics
353  // Data t0 from Calibrations_ftpc/ftpcElectronics
354 
355  TimeCoordinate += (tZero+deltat0)/micropertime;
356 
357  //LOG_DEBUG<<TimeCoordinate<<" "<<padpos<<" "<<(*x1)<<" "<<(*y1)<<endm;
358 
359  PadtransPerTimebin = (int) mParam->numberOfDriftSteps() / mDb->numberOfTimebins();
360  PadtransLower= (int) (TimeCoordinate*PadtransPerTimebin);
361 
362  float Rad=pradius[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
363  -(pradius[iRow + mDb->numberOfPadrowsPerSide()*(PadtransLower)]
364  -pradius[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)])/2;
365 
366 
367  PhiDeflect=pdeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
368  +(pdeflection[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)]
369  -pdeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower])/2;
370 
371  // Aenderungen damit Koord. richtig (vgl. chain+Frank)
372 
373  //PhiDeflect=0.0; // ? genauer vgl auch ZerofieldMaker !???
374 
375  float Phi;
376 
377  /* calculate phi angle from pad position */
378  // for FTPC West
379  if (iRow <10) {
380  Phi = mDb->radiansPerBoundary() / 2
381  + ((padpos-1) + 0.5) * mDb->radiansPerPad()
382  + PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
383  + mDb->radiansPerBoundary())+halfpi;
384  }
385 
386  //LOG_DEBUG<<"mDb->radiansPerBoundary() = "<<mDb->radiansPerBoundary()<<endm;
387  //LOG_DEBUG<<"mDb->radiansPerPad() = "<<mDb->radiansPerPad()<<endm;
388  //LOG_DEBUG<<"mDb->numberOfPads() = "<<mDb->numberOfPads()<<endm;
389 
390  // for FTPC East
391  /* Invert pad number (== Peak->PadPosition) for FTPC East */
392  /* (not yet understood where and why pad numbers were inverted) */
393  if (iRow>=10) {
394  Phi = mDb->radiansPerBoundary() / 2
395  + (159.5-(padpos-1)) * mDb->radiansPerPad()
396  -PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
397  + mDb->radiansPerBoundary())+halfpi;
398  }
399 
400  // debug
401  //LOG_DEBUG<<"====================================="<<endm;
402  //LOG_DEBUG<<"Radius = "<<Rad<<" | Phi = "<<Phi<<endl;
403  //LOG_DEBUG<<"====================================="<<endm;
404 
405  // Anm : Pointer genauer !!!!!!!??????????
406 
407  /* transform to cartesian */
408  (*x1) = Rad*cos(Phi);
409  if (iRow <10) {
410  (*x1) = -(*x1);
411  }
412 
413  (*y1) = Rad*sin(Phi);
414 
415  return kTRUE;
416 }
Generated by   doxygen 1.8.5