StRoot: Stl3Util/ftf/FtfSl3.cxx Source File

00001 /*:>-------------------------------------------------------------------
00002 **: FILE:     FtfSl3.cxx 
00003 **: HISTORY:
00004 **:           aug 23, 1999  add setup with number of hits and tracks
00005 **:           oct  6, 1999  ppy new version from Christof
00006 **:           oct 11, 1999  ppy call rawToLocal with variables instead of pointers
00007 **:           oct 11, 1999  ppy phiMin, phiMax changed to 0 and 2 pi
00008 **:           oct 21, 1999  ppy call to rawToGlobal with variables rather than pointers
00009 **:           oct 22, 1999  ppy back to pointers to make Christof happy
00010 **:
00011 **:           nov 24, 1999  cle extra method for filling US tracks deleted
00012 **:                         (FtfSl3::fillUSTracks) 
00013 **:                         added 'token' as argument for fill method
00014 **:                         changed formats to official L3 formats
00015 **:           dec  2, 1999  ppy setup and readSector modified
00016 **:                         sectorGeo class added to keep sector phase space
00017 **:                         readSector uses it to define tracker phase space
00018 **:           dec  6, 1999  ppy add xLastHit and yLastHit to type1_track
00019 **:           dec  6, 1999  ppy provisional patch to avoid generaring type1 tracks
00020 **:           dec  6, 1999  ppy pt in type 1 & 2 tracks is pt * charge  
00021 **:           dec  6, 1999  ppy method CanItBeMerged added
00022 **:           dec 16, 1999  ppy readSector and readMezzaninne look at
00023 **:                             byte order now
00024 **:           dec 20, 1999  ppy sectorGeo[23].phiMax set to correct value in setup 
00025 **:           dec 21, 1999  ppy fillTracks=0 deleted, leave default value (1)
00026 **:           dec 21, 1999  ppy bug corrected phiShift = 0 zero for non-border sectors
00027 **:           dec 21, 1999  ppy ptHelixFit set to zero in setParameters
00028 **:           dec 21, 1999  ppy maxChi2Primary =10 in setParameters
00029 **:           jan 27, 2000  ppy canItMerged beefed up
00030 **:           feb 10, 2000  ppy add xyError and zError
00031 **:           feb 10, 2000  ppy add xyError and zError
00032 **:           feb 14, 2000  ppy track length filled for tracks type 2 and 3
00033 **:           feb 17, 2000  ppy check ratiox and ratioy are no division by zero
00034 **:           feb 24, 2000  ppy hit id in readMezzaninne should be nHits+counter
00035 **:           mar 09, 2000  ppy trackLength changed to length
00036 **:           mar 09, 2000  ppy some changes to accomodate void pointers
00037 **:           mar 11, 2000  ppy fill hits
00038 **:           mar 29, 2000  ppy add errors to tracks in daq format
00039 **:           apr 06, 2000  cle change L3_CLUSTER to l3_cluster
00040 **:           apr 08, 2000  ppy add variables to cut cluster for high/low charge/time bin
00041 **:           apr 18, 2000  ppy include readout board and mezzannine in readMezzanine    
00042 **:                             arguments
00043 **:           apr 19, 2000  cs  add dEdx from Christof
00044 **:           apr 19, 2000  cs  discard one-apd cluster (flag == 1) in readMezzanine(...)
00045 **:           jun 28, 2000  ppy replace sl3CoordinateTransform with class
00046 **:                             St_l3_Coordinate_Transformer
00047 **:           jul 12, 2000  ppy get rid of mergeFlag in fillTracks (not used)
00048 **:           aug 21, 2000  ppy replace printf with ftfLog
00049 **:<------------------------------------------------------------------*/
00050 #include "Stl3Util/ftf/FtfSl3.h"
00051 
00052 #include "Stl3Util/base/L3swap.h"
00053 
00054 //******************************************************************
00055 //   Check whether tracs are mergable
00056 //******************************************************************
00057 int FtfSl3::canItBeMerged ( FtfTrack* tTrack ) {
00058    double rTpcMin =  50. ;
00059    double rTpcMax = 210. ;
00060 
00061    if ( nHits < 15 ) return 0 ;
00062 //
00063 //   Get circle parameters
00064 //
00065 // double rc = tTrack->pt / ( 2.9979e-3 * para.bField );
00066 // double xc = para.xVertex - tTrack->a2Xy / ( 2. * tTrack->a1Xy ) ;
00067 // double yc = para.yVertex - 1.   /  ( 2. * tTrack->a1Xy ) ;
00068    double x0 = tTrack->r0 * cos (tTrack->phi0);
00069    double y0 = tTrack->r0 * sin (tTrack->phi0);
00070    double rc = tTrack->pt / (bFactor * para.bField);
00071    double trackPhi0 = tTrack->psi + tTrack->q * 0.5 * M_PI / fabs ((double) tTrack->q);
00072    double xc = x0 - rc * cos (trackPhi0);
00073    double yc = y0 - rc * sin (trackPhi0);
00074 //   
00075 //   Calculate intersection with tracking area (sector or supersector) boundaries
00076 //
00077    double localPhi[2] ;
00078    localPhi[0] = para.phiMin ;
00079    localPhi[1] = para.phiMax ;
00080 
00081 // ftfLog ( " mergable track pt %f \n", tTrack->pt ) ;
00082    double sinPhi, cosPhi, tanPhi ;  
00083    double a, b, c, b2minus4ac ;
00084    double x1, y1, x2, y2, r1, r2 ;
00085    c = xc * xc + yc * yc - rc * rc ;
00086    for ( int i = 0 ; i < 2 ; i++ ) {
00087       sinPhi = sin(localPhi[i]); 
00088       cosPhi = cos(localPhi[i]); 
00089       tanPhi = tan(localPhi[i]); 
00090       a = 1. + tanPhi * tanPhi ;
00091       b = -2. * xc - 2. * yc * tanPhi ;
00092 
00093       b2minus4ac = b * b - 4. * a * c ;
00094       if ( b2minus4ac > 0 ) {
00095          double rootB2Minus4ac = ::sqrt(b2minus4ac);
00096 
00097          x1 = 0.5 * (-b + rootB2Minus4ac) / a ;
00098          y1 = x1 * tanPhi ;
00099          r1 = ::sqrt(x1*x1+y1*y1);
00100 
00101 //       ftfLog ( " x1 y1 %e %e cos sin %e %e\n", x1, y1, cosPhi, sinPhi ) ; 
00102          double ratiox = 0 ;
00103          if ( cosPhi != 0 ) ratiox = x1/cosPhi ;
00104          double ratioy = 0 ;
00105          if ( sinPhi != 0 ) ratiox = y1/sinPhi ;
00106 //       ftfLog ( " ratiox ratioy %e %e \n", ratiox, ratioy ) ;
00107          if ( ratiox >= 0 && ratioy >= 0 && 
00108               r1 > rTpcMin   &&  r1 < rTpcMax ) {
00109 //          ftfLog ( "phi %f boundary crossed !!!!\n", localPhi[i] ) ;
00110             return 1 ;
00111          }
00112 
00113          x2 = 0.5 * (-b - rootB2Minus4ac) / a ;
00114          y2 = x2 * tanPhi ;
00115          r2 = ::sqrt(x2*x2+y2*y2);
00116          ratiox = 0 ;
00117          if ( cosPhi != 0 ) ratiox = x2/cosPhi ;
00118          ratioy = 0 ;
00119          if ( sinPhi != 0 ) ratioy = y2/sinPhi ;
00120 //       ftfLog ( " ratiox ratioy %e %e \n", ratiox, ratioy ) ;
00121 //       ftfLog ( " x2 y2 %e %e cos sin %e %e\n", x2, y2, cosPhi, sinPhi ) ; 
00122          if ( ratiox >= 0 && ratioy >= 0 &&  
00123               r2 > rTpcMin   &&  r2 < rTpcMax ) {
00124 //          ftfLog ( "phi %f boundary crossed !!!!\n", localPhi[i] ) ;
00125             return 1 ;
00126          }
00127       }
00128 
00129    }
00130    //
00131    //  Check whether particle cross membrane
00132    //
00133    double zMembrane = 0 ;
00134    double sToMembrane = (zMembrane - tTrack->z0 ) / tTrack->tanl ;   
00135    double dangle      =  sToMembrane / rc ;    
00136    double angle       = dangle + trackPhi0  ; 
00137    double xMembrane   = xc + rc * cos(angle) ;
00138    double yMembrane   = yc + rc * sin(angle) ;
00139    double rMembrane   = ::sqrt(xMembrane*xMembrane+yMembrane*yMembrane);
00140    if ( rMembrane > rTpcMin   &&  rMembrane < rTpcMax ) {
00141 //    ftfLog ( "Membrane crossed at r %f !!!!\n", rMembrane ) ;
00142       return 1 ;
00143    }
00144 //
00145    return 0 ;
00146 }
00147 //******************************************************************
00148 //   fill Hits
00149 //******************************************************************
00150 int FtfSl3::fillHits ( unsigned int maxBytes, char* buff, unsigned int token ) {
00151 //
00152 //   Check enough space
00153 //
00154    if ( (nHits-1) * sizeof(l3_cluster) + sizeof(L3_SECCD) > maxBytes ) {
00155       ftfLog ( "FtfSl3::fillHits: buffer with %d bytes too small \n", maxBytes ) ;
00156       return 0 ;
00157    }
00158 //
00159 //   Loop over tracks
00160 //
00161    L3_SECCD* head = (L3_SECCD *)buff ;
00162 
00163    head->nrClusters_in_sector = nHits ;
00164    // bankHeader
00165    memcpy(head->bh.bank_type,CHAR_L3_SECCD,8);
00166    head->bh.bank_id    = sectorNr;
00167    head->bh.format_ver = DAQ_RAW_FORMAT_VERSION ;
00168    head->bh.byte_order = DAQ_RAW_FORMAT_ORDER ;
00169    head->bh.format_number = 0;
00170    head->bh.token      = token;
00171    head->bh.w9         = DAQ_RAW_FORMAT_WORD9;
00172    head->bh.crc        = 0; //don't know yet....
00173 
00174     
00175    l3_cluster* hitP    = (l3_cluster *)(head->cluster);
00176    FtfHit *nextHit ;
00177    int i ;
00178    int counter = 0 ;
00179 
00180 
00181    for ( i = 0 ; i < nTracks ; i++ ) {
00182       for  ( nextHit = (FtfHit *)(track[i].firstHit) ;
00183             nextHit != 0 ;
00184             nextHit = (FtfHit *)nextHit->nextTrackHit) {
00185          hitP->padrow  = nextHit->row ;
00186          hitP->pad     = nextHit->buffer1 ;
00187          hitP->time    = nextHit->buffer2 ;
00188          hitP->flags   = nextHit->flags   ;
00189          hitP->charge  = (short)nextHit->q  ;
00190          hitP->RB_MZ   = (short)nextHit->hardwareId  ;
00191          hitP->trackId = track[i].id ;
00192          if ( !(nextHit->track) ) {
00193             ftfLog ( "FtfSl3:fillHits: we got a problem track %d row hit %d \n",
00194                       track[i].id, hitP->padrow ) ;
00195          }
00196          hitP++ ;
00197          counter++ ;
00198       }
00199    }
00200 //
00201 //   Loop over hits
00202 //
00203    for ( i = 0 ; i < nHits ; i++ ) {
00204       nextHit = &(hit[i]) ;
00205       if ( nextHit->track ) continue ; 
00206       hitP->padrow  = nextHit->row ;
00207       hitP->pad     = nextHit->buffer1 ;
00208       hitP->time    = nextHit->buffer2 ;
00209       hitP->flags   = nextHit->flags   ;
00210       hitP->charge  = (short)nextHit->q       ;
00211       hitP->RB_MZ   = (short)nextHit->hardwareId  ;
00212       hitP->trackId = 0 ;
00213       hitP++ ;
00214       counter++ ;
00215    }
00216 
00217    if ( counter != nHits ) {
00218       ftfLog ( "FtfSl3:fillHits: Warning! Actual number of hits written %d total nHits %d \n", 
00219                                  counter, nHits );
00220    }
00221 
00222    head->bh.length = ((char *)hitP-buff) / 4 ;
00223 
00224    return ((char *)hitP-buff) ;
00225 }
00226 //******************************************************************
00227 //   fill tracks
00228 //******************************************************************
00229 int FtfSl3::fillTracks ( int maxBytes, char* buff, unsigned int token ) {
00230 //
00231 //   Loop over tracks
00232 //
00233     int counter = 0 ;
00234 
00235     counter = nTracks ;
00236 //
00237     // this is done to ease calcultion of offsets further down.
00238     unsigned int headSize;
00239 
00240     if(counter) 
00241       headSize = sizeof(struct bankHeader); 
00242     else 
00243       headSize = 0; 
00244 
00245 
00246     char *pointer = buff + sizeof(struct L3_SECTP) ;  
00247 
00248     char *endOfData = pointer + headSize
00249                        + counter * sizeof(struct local_track) ; // end
00250 
00251     int nBytes = (endOfData-buff) ; 
00252     if ( nBytes > maxBytes ) {
00253        ftfLog ( "FtfSl3::fillTracks: %d bytes needed, %d max, too short a buffer \n ",  nBytes, maxBytes ) ;
00254        return -1 ;
00255     }
00256 //
00257 //  Fill header
00258 //
00259 
00260     struct L3_SECTP *head = (struct L3_SECTP *)buff ;
00261 
00262     // bankHeader
00263     //  head->bh.bank_type = CHAR_L3_SECTP;
00264     memcpy(head->bh.bank_type,CHAR_L3_SECTP,8);
00265     head->bh.length     = sizeof(struct L3_SECTP) / 4 ;
00266     head->bh.bank_id    = sectorNr;
00267     head->bh.format_ver = DAQ_RAW_FORMAT_VERSION ;
00268     head->bh.byte_order = DAQ_RAW_FORMAT_ORDER ;
00269     head->bh.format_number = 1; // 1 means only pointing to local_tracks
00270     head->bh.token      = token;
00271     head->bh.w9         = DAQ_RAW_FORMAT_WORD9;
00272     head->bh.crc        = 0; //don't know yet....
00273 
00274     head->nHits    = (unsigned int)nHits ;
00275     head->nTracks  = (unsigned int)nTracks ;
00276     head->cpuTime  = (unsigned int)rint(cpuTime * 1000000);   // to get microsec
00277     head->realTime = (unsigned int)rint(realTime * 1000000) ; // same here
00278     head->xVert  = (int)rint(para.xVertex * 1000000);   // to get cm*10**-6
00279     head->yVert  = (int)rint(para.yVertex * 1000000); 
00280     head->zVert  = (int)rint(para.zVertex * 1000000); 
00281     head->para = 1 ;
00282     head->banks[0].off = (pointer-buff) / 4;
00283     head->banks[0].len = (endOfData-pointer)/ 4 ;
00284     head->banks[1].off = 0 ;
00285     head->banks[1].len = 0 ;
00286     head->banks[2].off = 0 ;
00287     head->banks[2].len = 0 ;
00288 
00289 // done with L3_SECTP
00290 
00291 
00292 //  L3_STKXD banks:
00293 
00294     struct L3_LTD* LTD = (struct L3_LTD*)pointer;
00295 
00296     struct local_track* uSTrack = 
00297         (struct local_track *)(pointer + headSize) ;
00298 
00299 
00300     if(headSize) {
00301         memcpy(LTD->bh.bank_type, CHAR_L3_LTD, 8);
00302         LTD->bh.length     = (sizeof(struct bankHeader) + 
00303                                       counter * sizeof(local_track))/4 ;
00304         LTD->bh.bank_id    = sectorNr ;
00305         LTD->bh.format_ver = DAQ_RAW_FORMAT_VERSION ;
00306         LTD->bh.byte_order = DAQ_RAW_FORMAT_ORDER ;
00307         LTD->bh.format_number = 0 ;
00308         LTD->bh.token      = token ;
00309         LTD->bh.w9         = DAQ_RAW_FORMAT_WORD9 ;
00310         LTD->bh.crc = 0; // for now!!!
00311     }   
00312 //
00313     for ( int i = 0 ; i < nTracks ; i++ ) {
00314        //
00315        uSTrack->id    = track[i].id ;
00316        if ( canItBeMerged ( &(track[i]) ) ) {
00317           uSTrack->id *= -1 ;
00318        }
00319        uSTrack->nHits = track[i].nHits ;
00320        uSTrack->ndedx = track[i].nDedx;
00321        uSTrack->innerMostRow = track[i].innerMostRow    ; 
00322        uSTrack->outerMostRow = track[i].outerMostRow    ; 
00323        uSTrack->xy_chisq = (unsigned short)rint(10 * track[i].chi2[0]
00324              /float(track[i].nHits));
00325        uSTrack->sz_chisq = (unsigned short)rint(10 * track[i].chi2[1]
00326              /float(track[i].nHits));
00327        uSTrack->dedx  = track[i].dedx  ; 
00328        uSTrack->pt    = track[i].pt * float(track[i].q * para.bFieldPolarity) ; 
00329        uSTrack->psi   = track[i].psi   ; 
00330        uSTrack->tanl  = track[i].tanl  ; 
00331        uSTrack->r0    = track[i].r0    ; 
00332        uSTrack->phi0  = track[i].phi0  ; 
00333        uSTrack->z0    = track[i].z0    ; 
00334        uSTrack->trackLength  = track[i].length  ;
00335        uSTrack->dpt          = (short)(32768. * track[i].dpt / track[i].pt ) ; 
00336        uSTrack->dpsi         = track[i].CompressOver1(track[i].dpsi,track[i].psi);
00337        uSTrack->dtanl        = track[i].CompressOver1(64.*track[i].dtanl,track[i].tanl);
00338        uSTrack->dz0          = (short)(1024. * track[i].dz0 ) ;
00339        uSTrack++ ;
00340     }
00341 //  delete array
00342 
00343 //
00344     return nBytes ;
00345 //
00346 }
00347 //******************************************************************
00348 //   get nr of tracks
00349 //******************************************************************
00350 int FtfSl3::getNrTracks (void) {
00351    return nTracks;
00352 }           
00353 //******************************************************************
00354 //   Initialize sl3 tracker
00355 //******************************************************************
00356 int FtfSl3::setup ( int maxHitsIn, int maxTracksIn ) {
00357 //
00358 //    Set parameters
00359 //
00360  // ftfLog("setting up tracker\n");
00361   para.setDefaults ( ) ;
00362 //
00363 //    Set Parameters
00364 //
00365   setParameters ( ) ;
00366 //
00367 //    Set sector geometry parameters
00368 //
00369   double toRad = acos(-1.)/180. ;
00370   sectorGeometry* sectorG = (sectorGeometry *)sectorGeo ;
00371   sectorG[ 0].phiMin =  45. * toRad ;
00372   sectorG[ 1].phiMin =  15. * toRad ;
00373   sectorG[ 2].phiMin = 345. * toRad ;
00374   sectorG[ 3].phiMin = 315. * toRad ;
00375   sectorG[ 4].phiMin = 285. * toRad ;
00376   sectorG[ 5].phiMin = 255. * toRad ;
00377   sectorG[ 6].phiMin = 225. * toRad ;
00378   sectorG[ 7].phiMin = 195. * toRad ;
00379   sectorG[ 8].phiMin = 165. * toRad ;
00380   sectorG[ 9].phiMin = 135. * toRad ;
00381   sectorG[10].phiMin = 105. * toRad ;
00382   sectorG[11].phiMin =  75. * toRad ;
00383   sectorG[12].phiMin = 105. * toRad ;
00384   sectorG[13].phiMin = 135. * toRad ;
00385   sectorG[14].phiMin = 165. * toRad ;
00386   sectorG[15].phiMin = 195. * toRad ;
00387   sectorG[16].phiMin = 225. * toRad ;
00388   sectorG[17].phiMin = 255. * toRad ;
00389   sectorG[18].phiMin = 285. * toRad ;
00390   sectorG[19].phiMin = 315. * toRad ;
00391   sectorG[20].phiMin = 345. * toRad ;
00392   sectorG[21].phiMin =  15. * toRad ;
00393   sectorG[22].phiMin =  45. * toRad ;
00394   sectorG[23].phiMin =  75. * toRad ;
00395 //
00396   sectorG[ 0].phiMax =  75. * toRad ;
00397   sectorG[ 1].phiMax =  45. * toRad ;
00398   sectorG[ 2].phiMax =  15. * toRad ;
00399   sectorG[ 3].phiMax = 345. * toRad ;
00400   sectorG[ 4].phiMax = 315. * toRad ;
00401   sectorG[ 5].phiMax = 285. * toRad ;
00402   sectorG[ 6].phiMax = 255. * toRad ;
00403   sectorG[ 7].phiMax = 225. * toRad ;
00404   sectorG[ 8].phiMax = 195. * toRad ;
00405   sectorG[ 9].phiMax = 165. * toRad ;
00406   sectorG[10].phiMax = 135. * toRad ;
00407   sectorG[11].phiMax = 105. * toRad ;
00408   sectorG[12].phiMax = 135. * toRad ;
00409   sectorG[13].phiMax = 165. * toRad ;
00410   sectorG[14].phiMax = 195. * toRad ;
00411   sectorG[15].phiMax = 225. * toRad ;
00412   sectorG[16].phiMax = 255. * toRad ;
00413   sectorG[17].phiMax = 285. * toRad ;
00414   sectorG[18].phiMax = 315. * toRad ;
00415   sectorG[19].phiMax = 345. * toRad ;
00416   sectorG[20].phiMax =  15. * toRad ;
00417   sectorG[21].phiMax =  45. * toRad ;
00418   sectorG[22].phiMax =  75. * toRad ;
00419   sectorG[23].phiMax = 105. * toRad ;
00420 
00421   double etaMin = 0.4 ;
00422   double etaMax = 2.3 ;
00423 
00424   for ( int sector = 0 ; sector < NSECTORS ; sector++ ) {
00425      sectorG[sector].phiShift = 0 ;
00426      if ( sector < 12 ) {
00427         sectorG[sector].etaMin = -1. * etaMin ;
00428         sectorG[sector].etaMax = etaMax ;
00429      }
00430      else {
00431         sectorG[sector].etaMin = -1. * etaMax ;
00432         sectorG[sector].etaMax = etaMin ;
00433      }
00434   }
00435   sectorG[ 2].phiShift = 16. * toRad ;
00436   sectorG[20].phiShift = 16. * toRad ; 
00437 
00438   // Reset tracker
00439   reset ( ) ;
00440 
00441   // Allocate memory 
00442   maxHits    = maxHitsIn ;
00443   maxTracks  = maxTracksIn ;
00444   hit        = new FtfHit[maxHits] ;
00445   track      = new FtfTrack[maxTracks] ;
00446 
00447   para.phiMin =  0.F * pi / 180.F ;
00448   para.phiMax =360.F * pi / 180.F  ;
00449   para.etaMin =-0.4F ;
00450   para.etaMax = 2.2F ;
00451 #ifdef TRDEBUG
00452   para.trackDebug = 24 ;
00453   para.debugLevel =  1 ;
00454 #endif
00455   return 0;
00456 }
00457 
00458 
00459 //******************************************************************
00460 //    Read cluster data from TPCMZCLD bank (using daqFormats.h)
00461 //
00462 //    returns number of hits from this mezzanine
00463 //
00464 //
00465 //     author: christof, derived from tonko's 'sl3'
00466 //     ppy:    modified to look at byte order      
00467 //******************************************************************
00468 //#define TRDEBUG
00469 int FtfSl3::readMezzanine (int sector,        int readOutBoard,
00470                            int mezzanninneNr, struct TPCMZCLD_local *mzcld) {
00471 
00472    unsigned int *tmp32;
00473    int row, rows;
00474    int len;
00475    int i, j;
00476    int counter;
00477 
00478    FtfHit *hitP = &hit[nHits];
00479 
00480    // Prepare transformation
00481    St_l3_xyz_Coordinate XYZ(0,0,0) ;
00482    St_l3_ptrs_Coordinate PTRS(0,0,0,0) ;
00483 
00484 
00485    counter = 0;
00486    short swapByte = 0 ;
00487    if     ( !checkByteOrder(mzcld->bh.byte_order) ) swapByte = 1 ;
00488 
00489    tmp32 = mzcld->padrowFiller ;
00490 
00491    rows = mzcld->rows;
00492    if ( swapByte ) rows = swap32(rows);
00493 
00494   // ftfLog("FtfSl3::readMezzaninne:  Found %d rows.\n",rows) ;
00495 
00496    for (i=0; i<rows; i++) {
00497       row = *tmp32++ ;
00498       len = *tmp32++ ;
00499       if ( swapByte ) {
00500          row = swap32(row);
00501          len = swap32(len);
00502       }
00503       //ftfLog("        Row# %d: clusters %d\n",row,len) ;
00504 
00505       for ( j=0; j<len; j++) {
00506          double fp, ft ;
00507          unsigned short pad ;
00508          unsigned short time ;
00509 
00510          struct xt {
00511             unsigned short x ;
00512             unsigned short t ;
00513          } *xt ;
00514          struct c {
00515             unsigned short flags ;
00516             unsigned short charge ;
00517          } *c ;
00518 
00519          xt = (struct xt *) tmp32++ ;
00520          c = (struct c *) tmp32++ ;
00521 //
00522 //   Discard 1 pad clusters
00523 // cs: discard one-pad cluster marked by the cluster finder as such (set first bit)
00524          //ftfLog ( "flags %d \n", c->flags ) ;
00525          if ( (c->flags & 1) == 1) continue;
00526          
00527 
00528 
00529          if ( !swapByte ) {
00530             pad  = xt->x  ;
00531             time = xt->t ;
00532          }
00533          else {
00534             pad  = swap16(xt->x)  ;
00535             time = swap16(xt->t)  ;
00536          }
00537          fp = (double)pad / 64. ;
00538          ft = (double)time/ 64. ;
00539 
00540 // Cut on timebins and clustercharge added by cle 02/21/00
00541          if ( ft < minTimeBin ) continue ;
00542          if ( ft > maxTimeBin ) continue ;
00543          if ( c->charge < minClusterCharge ) continue ;
00544          if ( c->charge > maxClusterCharge ) continue ;
00545 #ifdef MDEBUG
00546          ftfLog ( "FtfSl3:readMezzaninne row %d \n", row ) ;
00547          ftfLog ( "ft %f min max %d %d \n", ft, minTimeBin, maxTimeBin ) ;
00548          ftfLog ( "charge %d min max %d %d \n", c->charge, minClusterCharge, maxClusterCharge ) ;
00549 #endif
00550 //
00551          //      printf("sector:%i RB:%i pad: %f time: %f row: %i\n",sector, readOutBoard,fp,ft,row);
00552          PTRS.Setptrs((double)fp, (double)ft,(double)row, (double)sector) ;
00553 
00554 //       St_l3_xyz_Coordinate local;
00555 //       ((St_l3_Coordinate_Transformer *)coordinateTransformer)->raw_to_local(PTRS,local) ;
00556 //       ((St_l3_Coordinate_Transformer *)coordinateTransformer)->local_to_global(PTRS,local,XYZ) ;
00557          ((St_l3_Coordinate_Transformer *)coordinateTransformer)->raw_to_global(PTRS,XYZ) ;
00558 
00559          hitP->id  = nHits+counter ;
00560          hitP->row = row ;
00561          hitP->sector = sector ;
00562 
00563          hitP->x   = (float) XYZ.Getx();
00564          hitP->y   = (float) XYZ.Gety();
00565          hitP->z   = (float) XYZ.Getz();
00566 //       ftfLog ( "x %f y %f z %f \n", hitP->x, hitP->y, hitP->z ) ;
00567 
00568          hitP->dx  = xyError ;
00569          hitP->dy  = xyError ;
00570          hitP->dz  = zError ;
00571          hitP->buffer1 = pad ;
00572          hitP->buffer2 = time ;
00573          if(embedded)
00574              hitP->flags   = (c->flags | (1<<7) );
00575          else
00576              hitP->flags   = c->flags ;
00577          hitP->q       = c->charge ;
00578          hitP->hardwareId = readOutBoard * 16 + mezzanninneNr ;
00579 
00580          hitP++;
00581 
00582          counter++;
00583 
00584          if ( (nHits+counter)>maxHits ) {
00585             ftfLog ( "Error - FtfSl3:read: Hit array too small: counter %d maxHits %d \n",
00586                   counter, maxHits ) ;
00587             return -1;
00588          }
00589       }
00590    }
00591    return counter;
00592 }
00593 
00594 //******************************************************************
00595 //     Read Sector from buffer (using daqFormats.h)
00596 //
00597 //
00598 //     author: christof, derived from tonko's 'sl3'
00599 //     ppy:    modified to look at byte order      
00600 //******************************************************************
00601 int FtfSl3::readSector ( struct TPCSECLP *seclp ) {
00602 
00603     struct TPCRBCLP *rbclp ;
00604     struct TPCMZCLD_local *mzcld ;
00605     int iRb, kMz;
00606     int sector; 
00607     int nHitsOfMz;
00608     embedded = 0;
00609     
00610     // reset sector-hit counter
00611     nHits = 0;
00612     //
00613     // check byte order of SECLP bank
00614     // byte swapping is needed
00615     short swapByte = 0 ;
00616     if     ( !checkByteOrder(seclp->bh.byte_order) ) swapByte = 1 ;
00617     
00618     nHitsOfMz = 0;
00619     
00620     sector = (unsigned int)seclp->bh.bank_id ;
00621     if ( swapByte ) sector = swap32(sector) ;
00622     
00623     sectorNr = sector;
00624     sectorGeometry* sectorG = (sectorGeometry *)sectorGeo ;
00625     para.phiMin = sectorG[sector-1].phiMin ;
00626     para.phiMax = sectorG[sector-1].phiMax ;
00627     para.etaMin = sectorG[sector-1].etaMin ;
00628     para.etaMax = sectorG[sector-1].etaMax ;
00629 
00630     //   Check Sector 
00631     if ( sectorNr < 1 && sectorNr > 24 ) {
00632         ftfLog ( "Error - FtfSl3::readSector: Wrong sector %d!\n",sectorNr);
00633         return -1 ;
00634     }
00635     
00636     // run over receiver boards
00637     for (iRb=0; iRb<SB_RB_NUM; iRb++) {
00638         
00639         if (iRb==6) { 
00640             sector++ ;  // the other Sector
00641             if ( sector == 4 ) { 
00642                 para.phiMin = 0. ;
00643                 para.phiMax = 1.1 ; // 60 degress a bit more than a rad
00644                 para.phiShift = 0.8 ; // ~45 degrees
00645             }
00646             else if ( sector == 22 ) {
00647                 para.phiMin = 0. ;
00648                 para.phiMax = 1.1 ; // 60 degress a bit more than a rad
00649                 para.phiShift = 0.27 ; // ~15 degrees
00650             }
00651             else { 
00652                 if ( sectorG[sector-1].phiMin < para.phiMin ) 
00653                     para.phiMin = sectorG[sector-1].phiMin ;
00654                 if ( sectorG[sector-1].phiMax > para.phiMax ) 
00655                     para.phiMax = sectorG[sector-1].phiMax ;
00656                 para.phiShift = 0 ;
00657             }
00658             if ( sectorG[sector-1].etaMin < para.etaMin ) 
00659                 para.etaMin = sectorG[sector-1].etaMin ;
00660             if ( sectorG[sector-1].etaMax > para.etaMax ) 
00661                 para.etaMax = sectorG[sector-1].etaMax ;
00662         }
00663 
00664         if ( !(unsigned int)seclp->rb[iRb].off) continue ;
00665         
00666         int off = (unsigned int)seclp->rb[iRb].off;
00667         if ( swapByte ) off = swap32(off);
00668         rbclp = (struct TPCRBCLP *)((char *)seclp + off * 4) ;
00669         int swapByteMezzaninne = 0 ;
00670         if ( !checkByteOrder(rbclp->bh.byte_order) ) swapByteMezzaninne = 1;
00671         
00672         // run over the 3 mezzanines
00673         for (kMz=0; kMz<RB_MZ_NUM; kMz++) {
00674             
00675             if(rbclp->mz[kMz].off) {
00676 #ifdef TRDEBUG
00677                 if ( debugLevel > 1 ) {
00678                     ftfLog ( "FtfSl3::readSector:       MZCLD %d exists!\n", 
00679                              kMz+1) ;
00680                 }
00681 #endif
00682             }
00683             else continue ;
00684             
00685             off = rbclp->mz[kMz].off ;
00686             if ( swapByteMezzaninne ) off = swap32(off);
00687             //       ftfLog ( "mezz %d off len %d %d \n",
00688             //                kMz, rbclp->mz[kMz].off, rbclp->mz[kMz].len ) ;
00689             mzcld = (struct TPCMZCLD_local *) ((char *)rbclp + off*4) ;
00690             
00691             if (mzcld) {
00692                 nHitsOfMz = readMezzanine (sector, iRb, kMz, mzcld);
00693                 if (nHitsOfMz<0) {
00694                     ftfLog ( "FtfSl3:readSector: wrong reading mezzanine \n" ) ;
00695                     return -1;
00696                 }
00697                 nHits += nHitsOfMz;
00698             }
00699             
00700         }
00701     }
00702     //
00703     //   Since phiMin/Max etaMin/etaMax may have changed reset tracker
00704     //
00705     reset();
00706     //
00707     return 0;
00708 } 
00709 
00710 
00711 
00712 
00713 
00714 
00715 
00716 
00717 
00718 
00719 //******************************************************************
00720 //     Read Sector from buffer (using daqFormats.h)
00721 //
00722 //     function to read in 2 sectors from different event and merge them 
00723 //     to be able to embed one event into another one
00724 //     
00725 //     
00726 //******************************************************************
00727 int FtfSl3::readSector ( struct TPCSECLP *seclp1, struct TPCSECLP *seclp2 ) {
00728 
00729    struct TPCRBCLP *rbclp ;
00730    struct TPCMZCLD_local *mzcld ;
00731    int iRb, kMz;
00732    int sector,tmpsector; 
00733    int nHitsOfMz;
00734    struct TPCSECLP *seclp; // 'generic' sector pointer
00735    short swapByte ; 
00736 
00737    // reset sector-hit counter
00738    nHits = 0;
00739    //
00740    // check byte order of SECLP bank
00741    // byte swapping is needed
00742    short swapByte1 = 0 ;
00743    short swapByte2 = 0;
00744    if     ( !checkByteOrder(seclp1->bh.byte_order) ) swapByte1 = 1 ;
00745    if     ( !checkByteOrder(seclp2->bh.byte_order) ) swapByte2 = 1 ;
00746 
00747    //   nHitsOfMz = 0;
00748 
00749    sector = (unsigned int)seclp1->bh.bank_id ;
00750    if ( swapByte1 ) sector = swap32(sector) ;
00751    tmpsector = (unsigned int)seclp2->bh.bank_id ;
00752    if ( swapByte2 ) tmpsector = swap32(tmpsector) ;
00753    // check that both sectors are consistent
00754    if(tmpsector != sector){
00755        ftfLog("Error in embedding got different sectors to read in\n");
00756        exit(1);
00757    }
00758    //ftfLog("sector: %i\n",sector);
00759 
00760    sectorNr = sector;
00761    sectorGeometry* sectorG = (sectorGeometry *)sectorGeo ;
00762    para.phiMin = sectorG[sector-1].phiMin ;
00763    para.phiMax = sectorG[sector-1].phiMax ;
00764    para.etaMin = sectorG[sector-1].etaMin ;
00765    para.etaMax = sectorG[sector-1].etaMax ;
00766    //
00767    //   Check Sector 
00768    //
00769    if ( sectorNr < 0 && sectorNr > 24 ) {
00770       ftfLog ( "Error - FtfSl3::readSector: Wrong sector %d!\n",sectorNr);
00771       return -1 ;
00772    }
00773 
00774    seclp = seclp1;
00775    swapByte = swapByte1;
00776    embedded = 0;
00777 
00778    for(int embed=0;embed<2;embed++){
00779        // two  loops first set the pointers to event one then after that is 
00780        // done to event 2, at end of for loop
00781        
00782        
00783        // run over receiver boards
00784        for (iRb=0; iRb<SB_RB_NUM; iRb++) {
00785            
00786            if (iRb==6) { 
00787                sector++ ;       // the other Sector
00788                if ( sector == 4 ) { 
00789                    para.phiMin = 0. ;
00790                    para.phiMax = 1.1 ; // 60 degress a bit more than a rad
00791                    para.phiShift = 0.8 ; // ~45 degrees
00792                }
00793                else if ( sector == 22 ) {
00794                    para.phiMin = 0. ;
00795                    para.phiMax = 1.1 ; // 60 degress a bit more than a rad
00796                    para.phiShift = 0.27 ; // ~15 degrees
00797                }
00798                else { 
00799                    if ( sectorG[sector-1].phiMin < para.phiMin ) para.phiMin = sectorG[sector-1].phiMin ;
00800                    if ( sectorG[sector-1].phiMax > para.phiMax ) para.phiMax = sectorG[sector-1].phiMax ;
00801                    para.phiShift = 0 ;
00802                }
00803                if ( sectorG[sector-1].etaMin < para.etaMin ) para.etaMin = sectorG[sector-1].etaMin ;
00804                if ( sectorG[sector-1].etaMax > para.etaMax ) para.etaMax = sectorG[sector-1].etaMax ;
00805            }
00806            //ftfLog ( "sector %d rb %d phi min/max %f %f \n", sector, iRb, para.phiMin, para.phiMax ) ;
00807            
00808            
00809            if ( !(unsigned int)seclp->rb[iRb].off) continue ;
00810            
00811            int off = (unsigned int)seclp->rb[iRb].off;
00812            if ( swapByte ) off = swap32(off);
00813            rbclp = (struct TPCRBCLP *)((char *)seclp + off * 4) ;
00814            int swapByteMezzaninne = 0 ;
00815            if     ( !checkByteOrder(rbclp->bh.byte_order) ) swapByteMezzaninne = 1 ;
00816            
00817            // run over the 3 mezzanines
00818            for (kMz=0; kMz<RB_MZ_NUM; kMz++) {
00819                
00820                if(rbclp->mz[kMz].off) {
00821 #ifdef TRDEBUG
00822                    if ( debugLevel > 1 ) {
00823                        ftfLog ( "FtfSl3::readSector:       MZCLD %d exists!\n", kMz+1) ;
00824                    }
00825 #endif
00826                }
00827                else continue ;
00828                
00829                off = rbclp->mz[kMz].off ;
00830                if ( swapByteMezzaninne ) off = swap32(off);
00831                //       ftfLog ( "mezz %d off len %d %d \n",
00832                //                kMz, rbclp->mz[kMz].off, rbclp->mz[kMz].len ) ;
00833                mzcld = (struct TPCMZCLD_local *) ((char *)rbclp + off*4) ;
00834                
00835                if (mzcld) {
00836                    nHitsOfMz = readMezzanine (sector, iRb, kMz, mzcld);
00837                    if (nHitsOfMz<0) {
00838                        ftfLog ( "FtfSl3:readSector: wrong reading mezzanine \n" ) ;
00839                        return -1;
00840                    }
00841                    nHits += nHitsOfMz;
00842                }
00843                
00844            }
00845        }
00846 
00847        seclp = seclp2;
00848        swapByte = swapByte2;
00849        embedded = 1;
00850        // reset sector
00851        sector = sectorNr;
00852    }
00853 
00854    //
00855    //   Since phiMin/Max etaMin/etaMax may have changed reset tracker
00856    //
00857    reset();
00858    //
00859    return 0;
00860 } 
00861 
00862 
00863 
00864 
00865 
00866 
00867 
00868 //*******************************************************************
00869 //    Process
00870 //*******************************************************************
00871 int FtfSl3::processSector ( ){ 
00872    //
00873    //   Reset hit track assignment
00874    //
00875    for ( int h = 0 ; h < nHits ; h++ ) {
00876       hit[h].track = 0 ;
00877    }
00878    para.eventReset = 1 ;
00879    nTracks         = 0 ;
00880    process ( ) ;
00881    if (para.dEdx) dEdx();
00882    //
00883    if ( debugLevel > 0 ) 
00884       ftfLog ( " FtfSl3::process: tracks %i Time: real %f cpu %f\n", 
00885             nTracks, realTime, cpuTime ) ;
00886    //
00887    //
00888    return 1;
00889 }
00890 
00891 //********************************************************************
00892 //     Calculates deposited Energy
00893 //********************************************************************
00894 int FtfSl3::dEdx ( ) {
00895 
00896    for ( int i = 0 ; i<nTracks ; i++ ){
00897          if (track[i].nHits<para.minHitsForDedx) {
00898                track[i].nDedx = 0;
00899                track[i].dedx  = 0;
00900                continue;
00901          }
00902          fDedx->TruncatedMean(&track[i]);
00903    }
00904   return 0;
00905 }
00906 //***************************************************************
00907 //   Set Default parameters
00908 //***************************************************************
00909 int FtfSl3::setParameters ( ) {
00910 
00911    // FtfPara* para = &(para) ;
00912 
00913    xyError                = 0.3 ;
00914    zError                 = 1.0 ;
00915    para.hitChi2Cut        =  50.F  ;
00916    para.goodHitChi2       =  10.F  ;
00917    para.trackChi2Cut      =  30.F  ;
00918    para.maxChi2Primary    = 0.   ;
00919    para.segmentRowSearchRange = 2   ;
00920    para.trackRowSearchRange = 3    ;
00921    para.dphi              = 0.08F  ;
00922    para.deta              = 0.08F  ;
00923    para.dphiMerge        = 0.01F  ;
00924    para.detaMerge        = 0.02F  ;
00925    para.etaMinTrack      = -2.2F  ;
00926    para.etaMaxTrack      =  2.2F  ;
00927 
00928    para.dEdx             =  1     ;
00929    // get Errors = 1 makes FPE's
00930    para.getErrors        =  0     ;
00931    para.goBackwards      =  1     ;
00932    para.goodDistance     =  5.F   ;
00933    para.mergePrimaries   =  0     ;
00934    para.maxDistanceSegment = 50.F ;
00935    para.minHitsPerTrack  = 5      ;
00936    para.nHitsForSegment  = 2      ;
00937    para.nEta             = 40     ;
00938    para.nEtaTrack        = 40     ;
00939    para.nPhi             = 10     ;
00940    para.nPhiTrack        = 40     ;
00941    para.nPrimaryPasses   = 1      ;
00942    para.nSecondaryPasses = 0      ;
00943    para.xyErrorScale     = 1.0F   ;
00944    para.szErrorScale     = 1.0F   ;
00945    para.phiClosed        = 0      ;
00946 
00947 //   para.ptMinHelixFit    = 100.F  ;
00948     para.ptMinHelixFit    = 0.  ;
00949    para.rVertex          = 0.F    ;
00950    para.xVertex          = 0.F    ;
00951    para.yVertex          = 0.F    ;
00952    para.zVertex          = 0.F    ;
00953    para.dxVertex         = 0.005F ;
00954    para.dyVertex         = 0.005F ;
00955    para.phiVertex        = 0.F    ;
00956 
00957    para.zMax             = 205. ;
00958    return 0;
00959 }
00960 
00961 
00962 
00963 
00964