StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
gl3Event.cxx
1 //: FILE: gl3Event.cc
3 //: HISTORY:
4 //: 3 dec 1999 version 1.00
5 //: 8 apr 2000 include modifications form Clemens
6 //: 6 jul 2000 add St_l3_Coordinate_Transformer
7 //: 12 jul 2000 merge tracks using parameters at closest
8 //: approach (from sl3) and them extrapolate
9 //: them to inner most hit for consumers down the line
10 //: 13 aug 2000 replace trackMerging with maxSectorNForTrackMerging
11 //:Jens Berger 03 jul 2001 makeVertex added: calcs vertex of the event
12 //:<------------------------------------------------------------------
13 #include "gl3Event.h"
14 #include <rtsLog.h>
15 #include "gl3Histo.h"
16 #include <DAQ_READER/daqReader.h>
17 #include "FtfSl3.h"
18 
19 #include <DAQ_READER/daq_dta.h>
20 #include <DAQ_TPC/daq_tpc.h>
21 #include <DAQ_TPX/daq_tpx.h>
22 #include <DAQ_SC/daq_sc.h>
23 
24 tpc_t *pTPC=NULL;
25 
26 //
27 // evp should already contain the event that we want to read
28 //
29 //
30 
31 int gl3Event::readFromEvpReader(daqReader *rdr, float bField)
32 {
33  daq_dta *dd;
34 
35  LOG(DBG, "Reader from EVP Reader: evt=%d token=%d",rdr->seq,rdr->token);
36 
37  // Clear this event...
38  resetEvent();
39  nHits = 0;
40 
41  // Setupt the magnetic field
42  if(bField == 1000) { // try to read from file
43  bField = -.5; // if all else fails, set to reverse full field!
44 
45  dd = rdr->det("sc")->get("legacy");
46  if(dd) {
47  dd->iterate();
48  sc_t *sc = (sc_t *)dd->Void;
49  if(sc->valid) bField = sc->mag_field;
50  }
51  }
52  if(fabs(bField) < .1) bField = .1;
53  setBField(bField);
54  LOG(DBG, "bField set to %f",bField);
55 
56  // need a tracker...
57  //coordinateTransformer->Set_parameters_by_hand(0.581, 217.039, 201.138 ); // AuAu200
58  coordinateTransformer->Set_parameters_by_hand(0.581, 217.039 - 12.8, 201.138 - 12.8 );
59  //coordinateTransformer->Set_parameters_by_hand(0.6176, 200.668, 201.138 ); // 3.85GeV
60  //coordinateTransformer->LoadTPCLookupTable("/RTS/conf/L3/map.bin");
61 
62  FtfSl3 *tracker = new FtfSl3(coordinateTransformer, rdr);
63 
64  tracker->setup();
65  tracker->para.bField = fabs(bField);
66 
67  // +1 , -1 or 0 if zero field...
68  tracker->para.bFieldPolarity = (bField>0) ? 1 : -1;
69 
70 
72  tracker->setXyError(.12) ;
73  tracker->setZError(.24) ;
74  tracker->para.ptMinHelixFit = 0.;
75  tracker->para.maxChi2Primary = 0.;
76  tracker->para.trackChi2Cut = 10 ;
77  tracker->para.hitChi2Cut = 50 ;
78  tracker->para.goodHitChi2 = 20 ;
80 
81  tracker->reset();
82 
83  // need temporary track memory...
84  L3_SECTP *sectp = NULL;
85  sectp = (L3_SECTP *)malloc(szSECP_max);
86 
87  int i;
88  for(i=0;i<24;i++) {
89  if((i%2) == 0) {
90  tracker->nHits = 0;
91  tracker->setTrackingAngles(i+1); // only set angles by hypersector...
92  }
93  LOG(DBG, "READ TPC data for sector %d (0x%x)",i+1,rdr);
94 
95  // read in clusters...
96  LOG(DBG, "Reading clusters");
97  if(i != 100000) {
98  sectorFirstHit[i+1] = nHits;
99  readITPCClustersFromEvpReader(rdr, i+1);
100  readClustersFromEvpReader(rdr, i+1);
101  }
102 
103  // Do tracking...
104  LOG(DBG, "Tracking...");
105  tracker->setClustersFromGl3Event(this, i+1);
106  //tracker->readSectorFromEvpReader(i+1);
107 
108  int ret=0;
109  // only do tracking on full hypersectors...
110  if((i%2) == 1) {
111  ret = tracker->processSector();
112  //LOG("JEFF", "processSector returns %d", ret);
113  if(ret == 1) {
114  free(sectp);
115  delete tracker;
116  return 1;
117  }
118 
119 
120  ret = tracker->fillTracks(szSECP_max, (char *)sectp, 0);
122 
123  //LOG("JEFF", "SECP size = %d",sectp->bh.length*4 + sectp->banks[0].len*4);
124 
125  int n = readSectorTracks((char *)sectp);
126  //LOG("JEFF", "Got %d tracks",n);
127 
128  if(n < 0) {
129  LOG(WARN, "Error reading tracker: sector %d\n",i,0,0,0,0);
130  continue;
131  }
132  }
133  }
134 
135  // LOG(NOTE, "FINAL");
136  //LOG(INFO, "Finalizing...\n");y
137  finalizeReconstruction();
138  free(sectp);
139  // LOG(NOTE, "FINAL2");
140 
141 #ifdef EMC_LEGACY
142  // If calibrations not loaded nothing should happen here...
143  emc.readFromEvpReader(rdr, mem);
144 #endif
145 
146  delete tracker;
147  return 0;
148 }
149 
150 // Assume global tpc structure already filled....
151 // s = sector from 1
152 //
153 void gl3Event::readClustersFromEvpReader(daqReader *rdr, int sector)
154 {
155  daq_dta *dd;
156 
157  pTPC = NULL;
158 
159  // Read the data....
160  dd = rdr->det("tpx")->get("legacy",sector);
161  if(!dd) {
162  dd = rdr->det("tpc")->get("legacy", sector);
163  }
164 
165  if(dd) {
166  dd->iterate();
167  pTPC = (tpc_t *)dd->Void;
168  }
169  else {
170  LOG(DBG, "No data for sector %d check for TPC",sector);
171  return;
172  }
173 
174  LOG(DBG, "have clusters? %p %d",pTPC, pTPC->has_clusters);
175  if(!pTPC->has_clusters) return;
176 
177  for(int r=0;r<45;r++) {
178 
179  for(int j=0;j<pTPC->cl_counts[r];j++) {
180  tpc_cl *c = &pTPC->cl[r][j];
181 
182  gl3Hit *gl3c = &hit[nHits];
183  nHits++;
184 
185  l3_cluster sl3c;
186  sl3c.pad = (int)((c->p - 0.5) * 64);
187  sl3c.time = (int)((c->t - 0.5) * 64);
188  sl3c.charge = c->charge;
189  sl3c.flags = c->flags;
190  sl3c.padrow = r;
191  sl3c.RB_MZ = 0; // need to fake rb=0 so "set" doesn't change my sector
192  //c->p1;
193  //c->p2;
194  //c->t1;
195  //c->t2;
196 
197  gl3c->set(coordinateTransformer, sector, &sl3c);
198 
199  //printf("i=%d nHits=%d (%f %f %f) %f %f %f\n",i,nHits,
200  // sl3c.pad / 64.0,
201  // sl3c.time / 64.0,
202  // sl3c.charge,
203  // gl3c->getX(),gl3c->getY(),gl3c->getZ());
204 
205  }
206  }
207 }
208 
209  int gl3Event::readITPCClustersFromEvpReader(daqReader *rdr, int sector) {
210 
211  int ncld = 0;
212  daq_dta *dd;
213  dd = rdr->det("itpc")->get("cld", sector);
214  if(dd) {
215  while(dd->iterate()) {
216  ncld++;
217 
218  int padrow = dd->row;
219  int sec = dd->sec;
220 
221  for(u_int i=0;i<dd->ncontent;i++) {
222  float pad = dd->cld[i].pad;
223  float tb = dd->cld[i].tb;
224  int charge = dd->cld[i].charge;
225  int flags = dd->cld[i].flags;
226  int padrow = dd->row;
227 
228  gl3Hit *gl3c = &hit[nHits];
229  nHits++;
230 
231  gl3c->setITPCHit(coordinateTransformer, sec, padrow, pad, tb, charge, flags);
232  }
233  }
234  }
235 
236  return ncld;
237  }
238 
239 
240 //####################################################################
241 // Constructor and Destructor
242 //####################################################################
243 
244 gl3Event::gl3Event( l3CoordinateTransformer* inTrans,
245  l3EmcCalibration* inBemcCalib,
246  l3EmcCalibration* inEemcCalib,
247  int mxHits, int mxTracks)
248  : emc(inBemcCalib, inEemcCalib)
249 {
250  hit = NULL;
251  track = NULL;
252  busy = 0;
253 
254  trackContainer = 0;
255  trackIndex = 0;
256  hitProcessing = 0;
257  maxSectorNForTrackMerging = 1000000;
258  coordinateTransformer = inTrans;
259 
260  setup( mxHits, mxTracks );
261  resetEvent();
262 };
263 
264 
265 gl3Event::~gl3Event( )
266 {
267  if ( hit != 0 ) delete[] hit ;
268  if ( track != 0 ) delete[] track;
269  if ( trackContainer != 0 ) delete[] trackContainer;
270  if ( trackIndex != 0 ) delete[] trackIndex ;
271 };
272 
273 
274 
275 
276 //####################################################################
277 // Setters and getters
278 //####################################################################
279 gl3Track* gl3Event::getTrack ( int n ) {
280  if ( n < 0 || n > nTracks ) {
281  fprintf ( stderr, " %d track index out of range \n", n );
282  return NULL;
283  }
284  return &(track[n]);
285 }
286 
287 
288 gl3Hit* gl3Event::getHit ( int n ) {
289  if ( n < 0 || n > nHits ) {
290  fprintf ( stderr, " %d hit index out of range \n", n );
291  return NULL;
292  }
293  return &(hit[n]);
294 }
295 
296 gl3Sector* gl3Event::getSector ( int n ) {
297  if ( n < 0 || n > nSectors ) {
298  fprintf ( stderr, " %d sector index out of range \n", n );
299  return NULL;
300  }
301  return &(sectorInfo[n]);
302 }
303 
304 
305 int gl3Event::getTrgCmd()
306 {
307  //return trgData.trgCmd;
308  return -1;
309 };
310 
311 int gl3Event::getTrgWord()
312 {
313  return trgData.triggerWord;
314 };
315 
316 int gl3Event::getZDC(int n)
317 {
318  return trgData.ZDC[n];
319 };
320 
321 int gl3Event::getCTB(int n)
322 {
323  return trgData.CTB[n];
324 };
325 
326 double gl3Event::getZDCVertex()
327 {
328  return ((double)(trgData.ZDC[9] - trgData.ZDC[8]) + 21.3) * 3.3;
329 };
330 
331 
332 // Getters for the 64bit bunch crossing number are available as
333 // 32 and 64 bit verions
334 // PLEASE CHECK!!!
335 unsigned int gl3Event::getBXingLo()
336 {
337  return trgData.bunchXing_lo;
338 }
339 
340 unsigned int gl3Event::getBXingHi()
341 {
342  return trgData.bunchXing_hi;
343 }
344 
345 unsigned long long gl3Event::getBXing()
346 {
347  unsigned long long bx_hi_long = trgData.bunchXing_hi;
348  unsigned long long bx_lo_long = trgData.bunchXing_lo;
349 
350  return (bx_hi_long << 32) | bx_lo_long;
351 };
352 
353 
354 
355 //####################################################################
356 // addTracks: take a pointer to the local_tracks of a sector and
357 // (possibly) merge them with the cuurently known tracks
358 //####################################################################
359 void gl3Event::addTracks ( short sector, int nTrk, local_track* localTrack ) {
360 //
361  gl3Track* lTrack = &(track[nTracks]) ;
362  local_track *trk = localTrack ;
363  int indexStore = -1 ;
364 //
365  int idTrack ;
366  for ( int i = 0 ; i < nTrk ; i++ ) {
367  lTrack->set ( sector, trk ) ;
368  lTrack->id = sector * 1000 + abs(trk->id) ;
369  lTrack->para = &para ;
370  lTrack->sector = sector ;
371  idTrack = trk->id ;
372  //
373  // Check where to store relation orig track<->final track
374  // but only if requested
375  //
376  if ( hitProcessing ) {
377  indexStore = -1 ;
378  if ( abs(idTrack) < maxTracksSector )
379  indexStore = (sector-1)*maxTracksSector + abs(idTrack) ;
380  else {
381  LOG(ERR, " gl3Event::addTracks: max number of tracks per Sector reached %d reached", idTrack ,0,0,0,0) ;
382  }
383  }
384  //
385  // if id from buffer is negative track is mergable
386  //
387  gl3Track* fatterTrack = 0 ;
388  if ( maxSectorNForTrackMerging > nTrk && idTrack < 0 ) {
389 
390  fatterTrack = lTrack->merge ( trackContainer ) ;
391  if ( fatterTrack ) {
392  if ( hitProcessing && indexStore > 0 ) {
393  trackIndex[indexStore] =
394  ((char *)fatterTrack - (char *)track )/sizeof(gl3Track)+1;
395  }
396  trk++ ;
397  nMergedTracks++ ;
398  continue ;
399  }
400  nMergableTracks++ ;
401  }
402 
403  // Store new track in trackIndex array
404  if ( hitProcessing && indexStore > 0 )
405  trackIndex[indexStore] = nTracks + 1;
406 
407  // Increment counters
408  lTrack++ ;
409  nTracks++ ;
410  trk++ ;
411  if ( nTracks+1 >= maxTracks ) {
412  LOG(ERR," gl3Event::addTracks: max number of tracks %d reached, sector: %i nrSectorTracks: %i", maxTracks, sector, nTrk ,0,0) ;
413  nTracks-- ;
414  lTrack--;
415  break;
416  }
417  }
418 }
419 
420 
421 
422 //####################################################################
423 // fillTracks: fill the gl3Tracks into the L3_GTD
424 //####################################################################
425 int gl3Event::fillTracks ( int maxBytes, char* buffer, unsigned int token ){
426  //
427  // Check enough space
428  //
429  int nBytesNeeded = sizeof(L3_GTD) + (nTracks-1) * sizeof(global_track) ;
430  if ( nBytesNeeded > maxBytes ) {
431  LOG(ERR, " gl3Event::writeTracks: %d bytes needed less than max = %d \n",
432  nBytesNeeded, maxBytes ,0,0,0) ;
433  return 0 ;
434  }
435 
436  L3_GTD* head = (L3_GTD *)buffer ;
437 
438  head->nHits = nHits;
439  head->xVert = vertex.Getx();
440  head->yVert = vertex.Gety();
441  head->zVert = vertex.Getz();
442  // bankHeader
443  //
444  memcpy(head->bh.bank_type,CHAR_L3_GTD,8);
445  head->bh.bank_id = 1;
446  head->bh.format_ver = DAQ_RAW_FORMAT_VERSION ;
447  head->bh.byte_order = DAQ_RAW_FORMAT_ORDER ;
448  head->bh.format_number = 0;
449  head->bh.token = token;
450  head->bh.w9 = DAQ_RAW_FORMAT_WORD9;
451  head->bh.crc = 0; //don't know yet....
452  head->bh.length = (sizeof(struct L3_GTD)
453  + (nTracks-1) * sizeof(struct global_track))/4 ;
454 
455 
456  //
457  // Loop over tracks now
458  //
459  global_track* oTrack = (global_track *)head->track ;
460  int counter = 0 ;
461  for ( int i = 0 ; i < nTracks ; i++ ) {
462  if ( fabs(track[i].z0) > 205 ) {
463 
464  nBadTracks++ ;
465  continue ;
466  }
467  oTrack->id = track[i].id ;
468  oTrack->flag = track[i].flag ;
469  oTrack->innerMostRow = track[i].innerMostRow ;
470  oTrack->outerMostRow = track[i].outerMostRow ;
471  oTrack->nHits = track[i].nHits ;
472  oTrack->ndedx = track[i].nDedx ;
473  oTrack->q = track[i].q ;
474  oTrack->chi2[0] = track[i].chi2[0] ;
475  oTrack->chi2[1] = track[i].chi2[1] ;
476  oTrack->dedx = track[i].dedx ;
477  oTrack->pt = track[i].pt ;
478  oTrack->phi0 = track[i].phi0 ;
479  oTrack->r0 = track[i].r0 ;
480  oTrack->z0 = track[i].z0 ;
481  oTrack->psi = track[i].psi ;
482  oTrack->tanl = track[i].tanl ;
483  oTrack->length = track[i].length ;
484  oTrack->dpt = track[i].dpt ;
485  oTrack->dpsi = track[i].dpsi ;
486  oTrack->dz0 = track[i].dz0 ;
487  oTrack->dtanl = track[i].dtanl ;
488  oTrack++ ;
489  counter++ ;
490  }
491  head->nTracks = counter ;
492 //
493 // return number of bytes
494 //
495  return ((char *)oTrack-buffer) ;
496 }
497 
498 //####################################################################
499 // readL3Data: read the L3 contributions for an event. Currently
500 // includes TPC data, but SVT and FTPC data will also
501 // be in this data block at some time.
502 //####################################################################
503 int gl3Event::readL3Data( L3_P* header )
504 {
505 
506  char* buffer = (char *)header;
507  //L3_P* header = (L3_P *)buffer ;
508 
509  int length, offset ;
510  char* trackPointer ;
511  char* hitPointer ;
512 
513  resetEvent ( );
514  int i ;
515  L3_SECP* sectorP ;
516  for ( i = 0 ; i < nSectors ; i++ ) {
517  length = header->sector[i].len ;
518 
519  if ( length==0 ) continue ;
520 
521  offset = 4 * header->sector[i].off ;
522  sectorP = (L3_SECP *)&(buffer[offset]);
523 
524  trackPointer = (char *)sectorP + sectorP->trackp.off * 4 ;
525 
526 
527  int nSectorTracks = 0;
528  if (sectorP->trackp.off) {
529  nSectorTracks = readSectorTracks ( trackPointer ) ;
530 
531  if ( nSectorTracks < 0 ) {
532  LOG(ERR, "gl3Event:readEvent: error reading tracks, sector %d", i+1,0,0,0,0);
533  return -1 ;
534  }
535  }
536 
537  if ( hitProcessing && sectorP->sl3clusterp.off ) {
538  hitPointer = (char *)sectorP + sectorP->sl3clusterp.off * 4 ;
539  readSectorHits ( hitPointer, nSectorTracks ) ;
540  }
541 
542  }
543 
544  if(header->bh.format_number>=5 && header->trig.len){
545  //l3Log("Reading TRG data\n");
546  //trgData.read( (int*)header + header->trig.off );
547  trgData.readL3P(header);
548 
549  }
550 
551 
552  if( header->bh.format_number>=7 ){
553  //l3Log("Reading EMC data\n");
554  emc.readRawData(header);
555  } else {
556  //l3Log("Not reading EMC: format_number=%i, len=%i" ,
557  // header->bh.format_number, header->emc[0].len);
558  }
559 #ifdef EVENTDISPLAY
560  // For best merging (as least as 7/12/00) tracks
561  // are passed from sl3 to gl3 at point of closest approach
562  // the event viewer wants them(at least for now) at
563  // inner most point so we extrapolate to inner most hit radius
564 
565  double radius ;
566 
567  for ( int i = 0 ; i < nTracks ; i++ ) {
568 
569  radius = coordinateTransformer->
570  GetRadialDistanceAtRow(track[i].innerMostRow-1) ;
571 
572  track[i].updateToRadius ( radius ) ;
573 
574  // If track outside of TPC move radius since row
575  // radius is equal or larger than DistanceAtRow
576  if ( fabs(track[i].z0) > 205 ) track[i].updateToRadius ( radius+5. ) ;
577  if ( fabs(track[i].z0) > 205 ) {
578  LOG(ERR, "gl3Event:: problem after extrapolation id %d z0 %f",
579  track[i].id, track[i].z0 ,0,0,0) ;
580  }
581  }
582 #endif
583  //
584  // declare event as busy
585  //
586  busy = 1 ;
587 
588  return 0 ;
589 }
590 
591 
592 //####################################################################
593 // Do last reconstruction steps before analysis
594 //####################################################################
595 int gl3Event::finalizeReconstruction()
596 {
597 
598  // LOG(NOTE, "F1");
599  if (vertexFinder & 0x01) // internal calculation
600  makeVertex();
601 
602  //LOG(NOTE, "F2");
603 
604  if ((vertexFinder & 0x02) && lmv) { // low mult vertex finder
605  //LOG(NOTE, "F2.5");
606  lmv->makeVertex(this);
607  Ftf3DHit vtx = lmv->getVertex();
608 
609  lmVertex.Setx(vtx.x);
610  lmVertex.Sety(vtx.y);
611  lmVertex.Setz(vtx.z);
612  }
613 
614  //LOG(NOTE, "F3");
615  Ftf3DHit vertex_ftf;
616  vertex_ftf.x = vertex.Getx();
617  vertex_ftf.y = vertex.Gety();
618  vertex_ftf.z = vertex.Getz();
619 
620  // Compute DCAs for all tracks
621  for (int i=0 ; i<getNTracks() ; i++) {
622  getTrack(i)->setDca(vertex_ftf);
623  }
624 
625  //LOG(NOTE, "F4");
626  return 0;
627 }
628 
629 
630 
631 //####################################################################
632 // readSectorHits: does what you expect
633 //####################################################################
634 int gl3Event::readSectorHits ( char* buffer, int nSectorTracks ){
635  L3_SECCD* head = (L3_SECCD *)buffer ;
636  //
637  // Check coordinate transformer is there
638  //
639  if ( !coordinateTransformer ) {
640  LOG(ERR, "gl3Event::readSectorHits: there is not Coordinate Transformer",0,0,0,0,0);
641  return 0 ;
642  }
643 
644 
645  //
646  // Check bank header type
647  //
648  if ( strncmp(head->bh.bank_type,CHAR_L3_SECCD,8) ) {
649  LOG(ERR, "gl3Event::readSectorHits: wrong bank type %s",
650  head->bh.bank_type,0,0,0,0 ) ;
651  LOG(ERR, " correct bank type would be %s", CHAR_L3_SECCD,0,0,0,0 ) ;
652  return 0 ;
653  }
654  int sector = head->bh.bank_id;
655  int nSectorHits = head->nrClusters_in_sector ;
656 
657  //
658  // Check number of hits
659  //
660  if ( nHits + nSectorHits > maxHits ) {
661  LOG(ERR, "gl3Event:readSectorHits: not enough space for hits in sector %d", sector,0,0,0,0 ) ;
662  LOG(ERR, " maxHits %d nSectorHits %d nHits %d", maxHits,
663  nSectorHits, nHits ,0,0) ;
664  return 0 ;
665  }
666 
667  l3_cluster* cluster = (l3_cluster *)head->cluster ;
668  l3_cluster* hitP ;
669  gl3Hit* gHitP = 0 ;
670 
671  for ( int i = 0 ; i < nSectorHits ; i++ ) {
672  hitP = &(cluster[i]) ;
673  //
674  // Only if hits are going to be used for analysis unpack them
675  //
676  if ( hitProcessing > 1 ) {
677  gHitP = &(hit[nHits+i]);
678  gHitP->set (coordinateTransformer, sector, hitP);
679  }
680  //
681  // Now logic to reset trackId in hits
682  // This is to take care of track merging
683  //
684  int trkId = hitP->trackId ;
685  if ( trkId < 0 || trkId > nSectorTracks ) {
686  LOG(ERR, "gl3Event:readSectorHits: %d wrong track id in hit of sector %d \n",
687  trkId, sector ,0,0,0) ;
688  continue ;
689  }
690  int indexStore = (sector-1)*maxTracksSector+trkId ;
691  if ( indexStore < 0 || indexStore > nSectors*maxTracksSector ) {
692  LOG(ERR, "gl3Event:readSectorHits: %d wrong indexStore\n",
693  indexStore ,0,0,0,0) ;
694  continue ;
695  }
696  int index = trackIndex[indexStore] - 1 ;
697  if ( index < 0 || index > nTracks ) continue ;
698  //
699  // Only if hits are gonig to be used the
700  // linked lists are set
701  //
702  if ( hitProcessing > 1 ) {
703  if ( track[index].firstHit == 0 )
704  track[index].firstHit = (void *)gHitP ;
705  else
706  ((gl3Hit *)(track[index].lastHit))->nextHit = (void *)gHitP ;
707  track[index].lastHit = (void *)gHitP ;
708  gHitP->trackId = track[index].id ;
709  }
710  //
711  // Modify trackId of clusters got from sl3
712  //
713  hitP->trackId = track[index].id ;
714  // l3Log ( "hit trackId %d \n", track[index].id ) ;
715 
716  }
717  nHits += nSectorHits ;
718 
719  return nSectorHits ;
720 }
721 
722 
723 //####################################################################
724 // readSectorTracks: guess what it does ;)
725 // fills some general info and calls addTracks()
726 //####################################################################
727 int gl3Event::readSectorTracks ( char* buffer ){
728 
729  struct L3_SECTP *head = (struct L3_SECTP *)buffer ;
730 
731  if ( strncmp(head->bh.bank_type,CHAR_L3_SECTP,8) ) {
732  LOG(ERR, "gl3Event::readSectorTracks, wrong bank type %s\n",
733  head->bh.bank_type,0,0,0,0 ) ;
734  return -1 ;
735  }
736 
737  int sector = head->bh.bank_id ;
738  if ( sector < 0 || sector > nSectors ) {
739  LOG(ERR," gl3Event::readSector: %d wrong sector \n", sector ,0,0,0,0) ;
740  return 1 ;
741  }
742 
743  gl3Sector* sectorP = &(sectorInfo[sector-1]) ;
744  sectorP->filled = 1 ;
745  sectorP->nHits = head->nHits ;
746  sectorP->nTracks = head->nTracks ;
747  sectorP->cpuTime = head->cpuTime ;
748  sectorP->realTime = head->realTime ;
749  sectorP->xVert = float(head->xVert)/1000000 ;
750  sectorP->yVert = float(head->yVert)/1000000 ;
751  sectorP->zVert = float(head->zVert)/1000000 ;
752  sectorP->rVert = sqrt((double)( sectorP->xVert*sectorP->xVert +
753  sectorP->yVert*sectorP->yVert)) ;
754  sectorP->phiVert = atan2((double)sectorP->yVert,(double)sectorP->xVert) ;
755  if ( sectorP->phiVert < 0 ) sectorP->phiVert += 2. * M_PI ;
756 
757  //sectorP->print();
758  //
759  // Set vertex parameters for track merging
760  //
761  para.xVertex = sectorP->xVert ;
762  para.yVertex = sectorP->yVert ;
763  para.zVertex = sectorP->zVert ;
764  para.rVertex = sectorP->rVert ;
765  para.phiVertex = sectorP->phiVert ;
766  //
767  char* pointer = head->banks[0].off * 4 + buffer ;
768  int nSectorTracks ;
769  //
770  if ( (head->banks[0].len > 0) && (head->bh.format_number > 0) ) {
771  // l3Log ( "banks[0].len %d\n", head->banks[0].len ) ;
772  nSectorTracks = (4 * head->banks[0].len - sizeof(struct bankHeader))
773  /sizeof(struct local_track);
774  }
775  else nSectorTracks = 0 ;
776  //
777  // Add tracks
778  //
779  if ( nSectorTracks > 0 ) {
780  struct L3_LTD* headerLocal = (struct L3_LTD*)pointer ;
781  local_track* localTrack = headerLocal->track ;
782  addTracks ( sector, nSectorTracks, localTrack ) ;
783  }
784 
785  //
786  return sectorP->nTracks ;
787 }
788 
789 
790 //####################################################################
791 // Reconstrucht the vertex and store it in gl3Event::vertex
792 //####################################################################
793 int gl3Event::makeVertex (){
794  // debug
795  //printf ( "doing gl3Vertex process!!!\n");
796 
797  // init
798  //short sector = 0 ;
799  gl3Track* gTrack ;
800  Ftf3DHit closestHit ;
801 
802  hVz->Reset();
803  hVx->Reset();
804  hVy->Reset();
805 
806  // Comment to use the vertex of the last event
807  vertex.Setxyz(0.0,0.0,0.0);
808 
809 
810  for(int iter = 0 ; iter<2; iter++ ) {
811  // loop over gtracks
812  for(int trkcnt = 0 ; trkcnt<getNTracks(); trkcnt++ ) {
813  gTrack = getTrack(trkcnt);
814 
815  //printf("-----track %d\n", trkcnt);
816 
817  // acept only tracks with nHits more then minNoOfHitsOnTrack
818  if ( gTrack->nHits > minNoOfHitsOnTrackUsedForVertexCalc &&
819  gTrack->pt > minMomUsedForVertexCalc) {
820  // bad bad bad baby! Wouldn't it be nicer to use Vx and Vz!
821  closestHit = gTrack->closestApproach(getVertex().Getx(),
822  getVertex().Gety());
823 
824 
825  //printf("---------- (%4.2f %4.2f %4.2f)\n",closestHit.x,closestHit.y,closestHit.z);
826  // fill histos with the coordinates of the closest point to x=y=0
827  hVz->Fill(closestHit.z,1.0);
828  hVx->Fill(closestHit.x,1.0);
829  hVy->Fill(closestHit.y,1.0);
830  }
831  } // for(int trkcnt = 0 ; trkcnt<event->getNTracks(); trkcnt++ )
832 
833  // get and set the weighted mean
834  vertex.Setxyz(hVx->getWeightedMean(6.0),
835  hVy->getWeightedMean(6.0),
836  hVz->getWeightedMean(4.0));
837 
838  } //for(int iter = 0 ; iter<2; iter++ )
839 
840 
841  return 0;
842 }
843 
844 //####################################################################
845 //
846 //####################################################################
847 int gl3Event::resetEvent ( ){
848  nHits = 0 ;
849  nTracks = 0 ;
850  nMergedTracks = 0 ;
851  nMergableTracks = 0 ;
852  nBadTracks = 0 ;
853  busy = 0 ;
854 
855  memset(sectorFirstHit, 0, sizeof(sectorFirstHit));
856 
857  // Reset tracks
858  memset(trackContainer, 0,
859  para.nPhiTrackPlusOne*para.nEtaTrackPlusOne*sizeof(FtfContainer));
860 
861  // Reset hits
862  if ( hitProcessing ) {
863  memset ( trackIndex, 0, maxTracksSector*nSectors*sizeof(int) ) ;
864  delete[] hit;
865  hit = new gl3Hit[maxHits];
866  }
867 
868  // Reset trigger data
869  for (int i=0; i<16; i++)
870  trgData.ZDC[i] = 0;
871 
872  for (int i=0; i<240; i++)
873  trgData.CTB[i] = 0;
874 
875  /*
876  for ( int i = 0 ; i < nTracks ; i++ ) {
877  track[i].firstHit = 0 ;
878  track[i].lastHit = 0 ;
879  track[i].nextTrack = 0 ;
880  }
881  */
882 
883  // Reset vertex seed
884  vertex.Setxyz(0.0, 0.0, 0.0);
885 
886  emc.reset();
887 
888  return 0 ;
889 }
890 
891 //####################################################################
892 //
893 //####################################################################
894 int gl3Event::setup ( int mxHits, int mxTracks )
895 {
896 
897  if ( mxHits < 0 || mxHits > 1000000 ) {
898  LOG(ERR, " gl3Event::setup: maxHits %d out of range \n", maxHits,0,0,0,0 ) ;
899  mxHits = 500000 ;
900  }
901 
902  if ( mxTracks < 0 || mxTracks > 1000000 ) {
903  LOG(ERR, " gl3Event::setup: maxTracks %d out of range \n", maxTracks,0,0,0,0 );
904  mxTracks = 50000 ;
905  }
906 
907 
908  maxHits = mxHits ;
909  maxTracks = mxTracks ;
910  maxTracksSector = maxTracks*2/ nSectors ;
911  hit = new gl3Hit[maxHits] ;
912  track = new gl3Track[maxTracks] ;
913  trackIndex = new int[maxTracksSector*nSectors];
914  //
915  // Merging variables
916  //
917  nMergedTracks = 0 ;
918 
919  para.nPhiTrackPlusOne = para.nPhiTrack + 1 ;
920  para.nEtaTrackPlusOne = para.nEtaTrack + 1 ;
921  //-----------------------------------------------------------------------
922  // If needed calculate track area dimensions
923  //-----------------------------------------------------------------------
924  para.phiSliceTrack = (para.phiMaxTrack - para.phiMinTrack)/para.nPhiTrack;
925  para.etaSliceTrack = (para.etaMaxTrack - para.etaMinTrack)/para.nEtaTrack;
926 
927  int nTrackVolumes = para.nPhiTrackPlusOne* para.nEtaTrackPlusOne ;
928  trackContainer = new FtfContainer[nTrackVolumes];
929  if(trackContainer == NULL) {
930  LOG(ERR, "Problem with memory allocation... exiting\n",0,0,0,0,0) ;
931  return 1 ;
932  }
933  para.primaries = 1 ;
934  para.ptMinHelixFit = 1.e60 ;
935 
936  nTracks = 0 ;
937 
938  //-----------------------------------------------------------------------
939  // instanziate for vertex calc
940  //-----------------------------------------------------------------------
941  minNoOfHitsOnTrackUsedForVertexCalc=14; // points
942  minMomUsedForVertexCalc=0.25; // GeV
943 
944  char hid[50] ;
945  char title[100] ;
946 
947  strcpy ( hid, "Vertex_Vz" ) ;
948  strcpy ( title, "Vertex_Vz" ) ;
949  hVz = new gl3Histo ( hid, title, 400, -200., 200. ) ;
950 
951  strcpy ( hid, "Vertex_Vx" ) ;
952  strcpy ( title, "Vertex_Vx" ) ;
953  hVx = new gl3Histo ( hid, title, 100,-10,10);
954 
955  strcpy ( hid, "Vertex_Vy" ) ;
956  strcpy ( title, "Vertex_Vy" ) ;
957  hVy = new gl3Histo ( hid, title, 100,-10,10);
958 
959  // -----------------------------------------------------------------------
960 
961  return 0 ;
962 }
void set(short sectorIn, local_track *trk)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: gl3Track.h:62
Definition: FtfSl3.h:49
int Reset()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: gl3Histo.cxx:121
int readFromEvpReader(daqReader *rdr, float bField=1000)
Definition: gl3Event.cxx:31
Definition: gl3Hit.h:24
Definition: daq_sc.h:6
Definition: daq_tpc.h:22
Definition: daq_tpc.h:11
void addTracks(bool cuts=false)
Add tracks to the list of the rendered objects from current MuDst event.
Definition: EdMu.C:161
int Fill(double x, double weight=1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ...
Definition: gl3Histo.cxx:53