StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
l3CoordinateTransformer.cxx
1 //:>-----------------------------------------------------------------
2 //: FILE: l3CoordinateTransformer.cxx
3 //: HISTORY:
4 //: may2000 version 1.00
5 //: 29jun2000 ppy SectorSin and SectorCos values changed
6 //: 29jun2000 ppy local_to_global changed
7 //: 29jun2000 ppy global_to_raw(int sector, int row, ... ) added
8 //: 29jun2000 ppy local_to_raw (int row, ... ) added
9 //:<------------------------------------------------------------------
10 //:>------------------------------------------------------------------
11 //: CLASS: l3CoordinateTransformer
12 //: DESCRIPTION: Transforms coordinates from/to: global, local and raw
13 //: AUTHOR: dfl - Dominik Flierl, flierl@bnl.gov
14 //:>------------------------------------------------------------------
15 
16 #include <stdio.h>
17 #include <iostream>
18 #include <iomanip>
19 #include <stdlib.h>
20 #include <string.h>
21 #include "l3CoordinateTransformer.h"
22 
23 #include <unistd.h>
24 #include <sys/mman.h>
25 #include <errno.h>
26 
27 #include <sys/types.h>
28 #include <sys/stat.h>
29 #include <fcntl.h>
30 #include <rtsLog.h>
31 
32 using namespace std;
33 
34 
35 l3CoordinateTransformer::l3CoordinateTransformer()
36 {
37  // initialize transformations
38  // reset start values
39  //max_tb_inner =0;
40  //max_tb_outer =0;
41  transformation_errors =0;
42 
43  //Use_transformation_provided_by_db() ;
44  Set_parameters_by_hand() ;
45  //Get_parameters_from_db() ;
46  //Use_transformation_provided_by_db() ;
47  //Print_parameters() ;
48 
49 
50  TPCmap = NULL;
51  //LoadTPCLookupTable();
52 }
53 
54 
55 l3CoordinateTransformer::~l3CoordinateTransformer()
56 {
57  // if (transformation_errors>1000)
58  // {
59  // cout << transformation_errors
60  // << " transformation errors occured.\n";
61  // }
62 
63  if(TPCmap) free(TPCmap);
64  TPCmap = NULL;
65 }
66 
67 #include "mapbin.h"
68 
69 int l3CoordinateTransformer::LoadTPCLookupTable(char *)
70 {
71  // if (TPCmap) free(TPCmap);
72  // TPCmap = NULL;
73 
74  // // open and map file
75  // int fd = open(mapfilename, O_RDONLY);
76 
77  // if (fd == -1) {
78  // LOG(ERR, "Unable to open transformation map '%s'.\n",
79  // mapfilename);
80 
81  // return -1;
82  // }
83 
84  // int filesize = lseek(fd, 0, SEEK_END);
85  // void *file = mmap(0, filesize, PROT_READ, MAP_PRIVATE, fd, 0);
86  // if ((int)file == -1) {
87  // LOG(ERR,"Unable to mmap transformation map '%s'.Aborting.\n",
88  // mapfilename);
89  // return -1;
90  // }
91 
92  // alread initialized...
93  if(TPCmap) return 0;
94 
95  void *file = (void *)mapdotbinfile;
96 
97 
98  int type = ((int *)file)[0];
99  if ((type != 3) && (type != 100)) { //only support conversion map types
100  // 3 and 100 (local coordinates)
101  LOG(ERR, "No valid map found header file",0,0,0,0,0);
102  return -1;
103  }
104 
105  int headerSize = ((int *)file)[1];
106 
107  dpad = ((float *)file)[3];
108  dtb = ((float *)file)[4];
109  maxtb = ((float *)file)[5];
110 
111  npad = (int) ceil( 182. / dpad);
112  ntb = (int) ceil( maxtb / dtb);
113 
114  TPCmap = (float *)malloc( MAPDOTBINFILE_SZ - headerSize*4);
115  if (TPCmap == NULL) {
116  LOG(ERR, "Cannot allocate memory for lookup table.\n");
117  return -1;
118  } else if (TPCmap == (float*)-1) {
119  LOG(ERR, "TPCmap=-1. What's that??? 0x%x", TPCmap);
120  }
121 
122  memcpy(TPCmap, ((float *)file)+headerSize, MAPDOTBINFILE_SZ - headerSize*4);
123 
124  return 0;
125 }
126 
127 void l3CoordinateTransformer::raw_to_global(const l3ptrsCoordinate &raw,
128  l3xyzCoordinate &global)
129 {
130  l3xyzCoordinate local;
131 
132  raw_to_local(raw, local);
133  local_to_global(raw, local, global);
134 }
135 
136 
137 /*
138 //______________________________
139 void l3CoordinateTransformer::raw_to_global(const l3ptrsCoordinate &raw ,
140  l3xyzCoordinate &global)
141 {
142 
143  typedef float binmap_t[45][npad+1][ntb+1][3];
144  binmap_t *binmap = (binmap_t *)TPCmap;
145 
146 
147  // grid coordinates
148  int ipad = (int)floor(raw.Getp()/dpad);
149  int itb = (int)floor(raw.Gett()/dtb);
150 
151  // position in grid square
152  float wpad = raw.Getp()/dpad - (float)ipad;
153  float wtb = raw.Gett()/dtb - (float)itb;
154 
155  int sec = (int)raw.Gets() - 1;
156  int row = (int)raw.Getr() - 1;
157 
158  float x1 = binmap[sec][row][ipad ][itb ][0];
159  float y1 = binmap[sec][row][ipad ][itb ][1];
160  float z1 = binmap[sec][row][ipad ][itb ][2];
161  float x2 = binmap[sec][row][ipad+1][itb ][0];
162  float y2 = binmap[sec][row][ipad+1][itb ][1];
163  float z2 = binmap[sec][row][ipad+1][itb ][2];
164  float x3 = binmap[sec][row][ipad+1][itb+1][0];
165  float y3 = binmap[sec][row][ipad+1][itb+1][1];
166  float z3 = binmap[sec][row][ipad+1][itb+1][2];
167  float x4 = binmap[sec][row][ipad ][itb+1][0];
168  float y4 = binmap[sec][row][ipad ][itb+1][1];
169  float z4 = binmap[sec][row][ipad ][itb+1][2];
170 
171 
172 
173  float x = (1-wpad)*(1-wtb)*x1 + wpad*(1-wtb)*x2
174  + wpad*wtb*x3 + (1-wpad)*wtb*x4;
175  float y = (1-wpad)*(1-wtb)*y1 + wpad*(1-wtb)*y2
176  + wpad*wtb*y3 + (1-wpad)*wtb*y4;
177  float z = (1-wpad)*(1-wtb)*z1 + wpad*(1-wtb)*z2
178  + wpad*wtb*z3 + (1-wpad)*wtb*z4;
179 
180  global.Setxyz(x,y,z);
181 
182 }
183 
184 */
185 
186 
187 
188 //______________________________
189 void l3CoordinateTransformer::raw_to_local(const l3ptrsCoordinate &raw ,
190  l3xyzCoordinate &local )
191 {
192  double pitch = (raw.Getr()<14) ? innerSectorPadPitch : outerSectorPadPitch;
193  double pads2move = raw.Getp() - numberOfPadsAtRow[(int)raw.Getr()-1]/2;
194  local.Setx( pitch * (pads2move - 0.5) ) ;
195  local.Sety( radialDistanceAtRow[(int)(raw.Getr())-1] ) ;
196 
197  // Timebucket_to_z different for inner and outer sector
198  if( raw.Getr() <= 13 )
199  {
200  local.Setz(drift_length_inner - raw.Gett() * lengthPerTb) ;
201  if (raw.Gett()>max_tb_inner) {
202  transformation_errors++;
203  }
204  }
205  else
206  {
207  local.Setz(drift_length_outer - raw.Gett() * lengthPerTb) ;
208  if (raw.Gett()>max_tb_outer) {
209  transformation_errors++;
210  }
211  }
212 }
213 //______________________________
214 void l3CoordinateTransformer::local_to_global(const l3ptrsCoordinate &raw ,
215  const l3xyzCoordinate &local,
216  l3xyzCoordinate &global )
217 {
218  int sector = (int) raw.Gets();
219 
220  // rotate local x,y coordinates
221  // 2x2 rotation matrix:
222  // ( cos b sin b )
223  // (-sin b cos b )
224  //
225  // caution: sector>12 needs x->-x and y->y (east side!)
226  double x = SectorCos[sector-1] * local.Getx() + SectorSin[sector-1] * local.Gety();
227  if (sector>12) x = -x;
228  double y = -1.* SectorSin[sector-1] * local.Getx() +
229  SectorCos[sector-1] * local.Gety();
230  double z = (sector<13) ? local.Getz() : -local.Getz() ;
231 
232  global.Setxyz(x,y,z);
233 
234 }
235 //______________________________
236 void l3CoordinateTransformer::global_to_raw(const l3xyzCoordinate &global,
237  l3ptrsCoordinate &raw )
238 {
239  l3xyzCoordinate local(0,0,0) ;
240  global_to_local( global, local, raw ) ;
241 
242 
243  local_to_raw( global ,local ,raw ) ;
244 }
245 //______________________________
246 void l3CoordinateTransformer::global_to_raw(int sector, int row,
247  const l3xyzCoordinate &global,
248  l3ptrsCoordinate &raw )
249 {
250  l3xyzCoordinate local(0,0,0) ;
251  global_to_local( sector, row, global, local ) ;
252  local_to_raw( row ,local ,raw ) ;
253 }
254 //______________________________
255 void l3CoordinateTransformer::global_to_local(const l3xyzCoordinate &global,
256  l3xyzCoordinate &local,
257  l3ptrsCoordinate &raw)
258 {
259 
260  // Get xyz right
261 
262  double y = global.Gety() ;
263  if (y== 0) y= 0.000000001;
264 
265  double x = 0;
266  if(global.Getz()>=0)
267  {
268  x = global.Getx() ;
269  local.Setz(global.Getz());
270  }
271  else
272  {
273  x = -(global.Getx()) ; // ATTENTION must be mirrowed for sectors 13 - 24 !!!
274  local.Setz(-(global.Getz()));
275  }
276 
277  // Prepare turn operation
278  double pi = 3.14159265358979323846;
279  double sec_border = tan(pi/12) ; // 15 degree
280  double turn_angle = -pi/6 ; // 30 degree
281  double sin_turn_angle = sin(turn_angle);
282  double cos_turn_angle = cos(turn_angle);
283  double sector = 0 ;
284 
285  if (y>=0 && fabs(x/y)<=sec_border)
286  {
287  // We are already in sector 12
288  sector = 12 ;
289  }
290  else
291  {
292  // We have to turn system until we are in first sector
293  while( y<0 || (fabs(x/y)>sec_border))
294  {
295  double xn = x*cos_turn_angle + y*sin_turn_angle ;
296  double yn = -x*sin_turn_angle + y*cos_turn_angle ;
297  x = xn ;
298  y = yn ;
299  sector++;
300  }
301  }
302 
303  // Set it
304  local.Setx(x);
305  local.Sety(y);
306 
307  if (global.Getz()<0)
308  {
309  raw.Sets(sector+12);
310  }
311  else
312  {
313  raw.Sets(sector);
314  }
315 
316 }
317 //______________________________
318 void l3CoordinateTransformer::global_to_local( int sector, int row,
319  const l3xyzCoordinate &global,
320  l3xyzCoordinate &local )
321 {
322 
323  // rotate global x,y coordinates back to local
324  // 2x2 rotation matrix:
325  // ( cos b -sin b )
326  // ( sin b cos b )
327  //
328  // caution: sector>12 needs x->-x and y->y (east side!)
329  double xGlobal = global.Getx();
330  if (sector>12) xGlobal = -xGlobal;
331 
332  double x = SectorCos[sector-1] * xGlobal - SectorSin[sector-1] * global.Gety() ;
333  double y = SectorSin[sector-1] * xGlobal + SectorCos[sector-1] * global.Gety() ;
334  double z = fabs(global.Getz());
335 
336  local.Setxyz(x,y,z);
337 
338  return ;
339 }
340 //______________________________
341 void l3CoordinateTransformer::local_to_raw(const l3xyzCoordinate &global,
342  const l3xyzCoordinate &local,
343  l3ptrsCoordinate &raw )
344 {
345  // first lets find the row
346  double y = local.Gety() ;
347  int row = 0;
348  int row_index = 0 ;
349  while( (fabs(radialDistanceAtRow[row_index]-y) > 0.5) && (row_index < 46) )
350  {
351  row_index++;
352  }
353  if (row_index==45 || y<59)
354  {
355  // no matching row found
356  //cerr << "Alert row not found !" << endl;
357  return ;
358  }
359  else
360  {
361  // yes row found !
362  row = row_index+1;
363  }
364 
365  // then lets go for the pad
366  double x = local.Getx();
367  double pitch = (row<=13) ? innerSectorPadPitch : outerSectorPadPitch ;
368  int half_num_pads_this_row = numberOfPadsAtRow[row-1]/2 ;
369  double pad = half_num_pads_this_row + x/pitch + 0.5 ;
370 
371  // finally lets get the bucket
372  double bucket = 0;
373  double z = local.Getz();
374  if (row<=13)
375  {
376  bucket = ( drift_length_inner - z )/lengthPerTb ;
377  if (z>drift_length_inner) {
378  transformation_errors++;
379  }
380  }
381  else
382  {
383  bucket = ( drift_length_outer - z )/lengthPerTb ;
384  if (z>drift_length_outer) {
385  transformation_errors++;
386  }
387  }
388 
389  // fill it
390  raw.Setp(pad);
391  raw.Sett(bucket);
392  raw.Setr(row);
393 }
394 //______________________________
395 void l3CoordinateTransformer::local_to_raw( int row ,
396  const l3xyzCoordinate &local ,
397  l3ptrsCoordinate &raw )
398 {
399  // then lets go for the pad
400  double x = local.Getx();
401  double pitch = (row<=13) ? innerSectorPadPitch : outerSectorPadPitch ;
402  int half_num_pads_this_row = numberOfPadsAtRow[row-1]/2 ;
403  double pad = half_num_pads_this_row + x/pitch + 0.5 ;
404 
405  // finally lets get the bucket
406  double bucket = 0;
407  double z = local.Getz();
408  if (row<=13)
409  {
410  bucket = ( drift_length_inner - z )/lengthPerTb ;
411  if (z>drift_length_inner) {
412  transformation_errors++;
413  }
414  }
415  else
416  {
417  bucket = ( drift_length_outer - z )/lengthPerTb ;
418  if (z>drift_length_outer) {
419  transformation_errors++;
420  }
421  }
422 
423  // fill it
424  raw.Setp(pad);
425  raw.Sett(bucket);
426  raw.Setr(row);
427 }
428 //______________________________
429 void l3CoordinateTransformer::Set_parameters_by_hand( const double mlengthPerTb,
430  const double mdrift_length_inner,
431  const double mdrift_length_outer)
432 {
433  lengthPerTb = mlengthPerTb ;
434  drift_length_inner = mdrift_length_inner ;
435  drift_length_outer = mdrift_length_outer ;
436 
437  // set max timebucket
438  max_tb_inner = drift_length_inner/lengthPerTb;
439  max_tb_outer = drift_length_outer/lengthPerTb;
440 }
441 
442 
443 //______________________________
444 void l3CoordinateTransformer::Use_transformation_provided_by_db()
445 {
446 #ifdef L3OFFLINE
447  // perform official transform and recalulate parameters
448  // the official transform doesn't take doubles for pad and time
449  // that's why we must extract the parameters ...
450  StTpcPadCoordinate padco;
451  StGlobalCoordinate glo;
452  StTpcCoordinateTransform tra(gStTpcDb);
453 
454  // shaping time
455  double tau3 = 3 * (gStTpcDb->Electronics()->tau() * 1e-09)
456  * (gStTpcDb->DriftVelocity()) ;
457 
458  // inner row (5) in sector 1
459  padco.setSector(1);
460  padco.setRow(5);
461  padco.setPad(1);
462  padco.setTimeBucket(100);
463 
464  tra(padco,glo);
465  double z_100 = glo.position().z();
466  padco.setTimeBucket(0);
467  tra(padco,glo);
468  double z_0 = glo.position().z();
469  double lengthPerTb_inner = fabs((z_0-z_100)/100) ;
470  drift_length_inner = fabs(glo.position().z()+tau3);
471 
472  //cout << "sector : " << padco.sector() << "\t";
473  //cout << "length per tb inner = " << lengthPerTb_inner << "\t";
474  //cout << "innerdrift : " <<drift_length_inner << endl;
475 
476  // outer row (30) in sector 1
477  padco.setSector(1);
478  padco.setRow(30);
479  padco.setPad(1);
480  padco.setTimeBucket(100);
481  tra(padco,glo);
482  z_100 = glo.position().z();
483  padco.setTimeBucket(0);
484  tra(padco,glo);
485  z_0 = glo.position().z();
486  double lengthPerTb_outer = fabs((z_0-z_100)/100) ;
487  drift_length_outer = fabs(glo.position().z()+tau3) ;
488 
489  //cout << "sector : " << padco.sector() << "\t";
490  //cout << "length per tb outer = " << lengthPerTb_outer << "\t";
491  //cout << "outerdrift : " << drift_length_outer << endl;
492 
493 
494  if (lengthPerTb_outer == lengthPerTb_inner )
495  {
496  lengthPerTb = lengthPerTb_outer ;
497  }
498  else
499  {
500  cerr << "Different driftvelocities in TPC observed -> Set to 0. " << endl;
501  lengthPerTb = 0;
502  }
503  //cout << "Constants set by using official transformation." << endl;
504 
505  // set max timebucket
506  max_tb_inner = drift_length_inner/lengthPerTb;
507  max_tb_outer = drift_length_outer/lengthPerTb;
508 
509 #else
510  cout << "This is not functional online.\n" ;
511 #endif
512 
513 };
514 
515 
516 //______________________________
517 void l3CoordinateTransformer::Get_parameters_from_db()
518 {
519 #ifdef L3OFFLINE
520  // connect to database
521  double driftvelocity = gStTpcDb->DriftVelocity();
522  double inner_effective_driftlength =
523  gStTpcDb->Dimensions()->innerEffectiveDriftDistance();
524  double outer_effective_driftlength =
525  gStTpcDb->Dimensions()->outerEffectiveDriftDistance();
526  double gatingrid = gStTpcDb->Dimensions()->gatingGridZ();
527  double t0pad = gStTpcDb->T0(1)->getT0(1,1);
528  double zinneroffset = gStTpcDb->Dimensions()->zInnerOffset();
529  double zouteroffset = gStTpcDb->Dimensions()->zOuterOffset();
530  double frequency = gStTpcDb->Electronics()->samplingFrequency();
531  double tzero = gStTpcDb->Electronics()->tZero();
532  double tau = gStTpcDb->Electronics()->tau();
533  double shapingtime = gStTpcDb->Electronics()->shapingTime();
534 
535  if (0)
536  {
537  cout << "Fresh from the db we got : \n" ;
538  cout << "The driftvelocity : " << driftvelocity << endl;
539  cout << "The effective innerdrift (not used) : "
540  << inner_effective_driftlength << endl ;
541  cout << "The effective outerdrift (not used) : "
542  << outer_effective_driftlength << endl ;
543  cout << "Innerzoffset : " << zinneroffset << endl ;
544  cout << "Outerzoffset : " << zouteroffset << endl ;
545  cout << "Gatinggrid : " << gatingrid << endl ;
546  cout << "t0pad : " << t0pad << endl ;
547  cout << "The tzero : " << tzero << endl;
548  cout << "Frequency ist set to : " << frequency << endl;
549  cout << "tau : " << tau << endl ;
550  cout << "shapingtime (not used) : " << shapingtime << endl ;
551  }
552 
553  // length per timebucket = 1 / frequency * driftvelocity
554  lengthPerTb = 1 / frequency / (1E6) * driftvelocity ; // cm per timebucket
555 
556  // driftlength = gatingrid - dv * t0 (this is a general t0) + dv/frequency * 0.5 (half a timebucket to get in the middle of a tb)
557  // + zoffset (differnt offset of inner and outer padrows) + t0 (per individual pad = 0 at 05/05/2000)
558  // + 3*tau*driftvelocity ( = shaping time (this depends on the used algorithm used ) this holds for weigthed mean see STARNOT
559  //
560  drift_length_inner = gatingrid-driftvelocity*1e-6*(-tzero+0.5/frequency)+
561  zinneroffset+t0pad+3*tau*1e-09*driftvelocity;
562  drift_length_outer = gatingrid-driftvelocity*1e-6*(-tzero+0.5/frequency)+
563  zouteroffset+t0pad+3*tau*1e-09*driftvelocity;
564 
565 
566  if (0)
567  {
568  cout << "inner drift length : " << drift_length_inner << endl;
569  cout << "outer drift length : " << drift_length_outer << endl;
570  cout << "lengthPerTb : " << lengthPerTb << endl;
571  }
572 
573  // set max timebucket
574  max_tb_inner = drift_length_inner/lengthPerTb;
575  max_tb_outer = drift_length_outer/lengthPerTb;
576 
577 #else
578  cout << "This is not functional online.\n" ;
579 #endif
580 
581 };
582 
583 
584 //______________________________
585 void l3CoordinateTransformer::Print_parameters()
586 {
587  cout << "l3CoordinateTransformer: Parameters used <== " << endl ;
588  cout << " Used parameters : " << endl ;
589  cout << " Length per tb : " << lengthPerTb << endl ;
590  cout << " drift_length_inner: " << drift_length_inner << endl ;
591  cout << " drift_length_outer: " << drift_length_outer << endl ;
592  cout << " max_tb_inner : " << max_tb_inner << endl ;
593  cout << " max_tb_outer : " << max_tb_outer << endl ;
594 
595 };
596 
598 // init our statics
600 double l3CoordinateTransformer::outerSectorPadPitch = 0.67; // cm
601 double l3CoordinateTransformer::innerSectorPadPitch = 0.335; // cm
602 
603 // number of pads in padrow
604 int l3CoordinateTransformer::numberOfPadsAtRow[] = {
605  88,96,104,112,118,126,134,142,150,158,166,174,182,
606  98,100,102,104,106,106,108,110,112,112,114,116,
607  118,120,122,122,124,126,128,128,130,132,134,136,
608  138,138,140,142,144,144,144,144
609 };
610 // radial distance (center pad) from detector center in cm
611 double l3CoordinateTransformer::radialDistanceAtRow[] = {
612  60.0, 64.8, 69.6, 74.4, 79.2, 84.0, 88.8, 93.60, // 7 * 4.80 cm spacing
613  98.8, 104., 109.20, 114.4, 119.6, // 5 * 5.20 cm spacing
614  127.195, 129.195, 131.195, 133.195, 135.195, // 32 * 2.00 cm spacing
615  137.195, 139.195, 141.195, 143.195, 145.195,
616  147.195, 149.195, 151.195, 153.195, 155.195,
617  157.195, 159.195, 161.195, 163.195, 165.195,
618  167.195, 169.195, 171.195, 173.195, 175.195,
619  177.195, 179.195, 181.195, 183.195, 185.195,
620  187.195, 189.195
621 };
622 // sector-rotation factors: 30 degree steps
623 double l3CoordinateTransformer::SectorSin[] = {
624  0.5, 0.866025404, // 1-2
625  1.0, 0.866025404, // 3-4
626  0.5, 0., // 5-6
627  -0.5, -0.866025404, // 7-8
628  -1.0, -0.866025404, // 9-10
629  -0.5, 0., // 11-12
630  0.5, 0.866025404, // 13-14
631  1.0, 0.866025404, // 15-16
632  0.5, 0., // 17-18
633  -0.5, -0.866025404, // 19-20
634  -1.0, -0.866025404, // 21-22
635  -0.5, 0. // 23-24
636 };
637 // sector-rotation factors: 30 degree steps
638 double l3CoordinateTransformer::SectorCos[] = {
639  0.866025404, 0.5, // 1-2
640  0., -0.5, // 3-4
641  -0.866025404, -1.0, // 5-6
642  -0.866025404, -0.5, // 7-8
643  0., 0.5, // 9-10
644  0.866025404, 1.0, // 11-12
645  0.866025404, 0.5, // 13-14
646  0., -0.5, // 15-16
647  -0.866025404, -1.0, // 17-18
648  -0.866025404, -0.5, // 19-20
649  0., 0.5, // 21-22
650  0.866025404, 1.0 // 23-24
651 };
652 
Generated by   doxygen 1.8.5