16 #include <DAQ_READER/daqReader.h>
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>
35 LOG(DBG,
"Reader from EVP Reader: evt=%d token=%d",rdr->seq,rdr->token);
45 dd = rdr->det(
"sc")->get(
"legacy");
49 if(sc->valid) bField = sc->mag_field;
52 if(fabs(bField) < .1) bField = .1;
54 LOG(DBG,
"bField set to %f",bField);
58 coordinateTransformer->Set_parameters_by_hand(0.581, 217.039 - 12.8, 201.138 - 12.8 );
62 FtfSl3 *tracker =
new FtfSl3(coordinateTransformer, rdr);
65 tracker->para.bField = fabs(bField);
68 tracker->para.bFieldPolarity = (bField>0) ? 1 : -1;
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 ;
85 sectp = (
L3_SECTP *)malloc(szSECP_max);
91 tracker->setTrackingAngles(i+1);
93 LOG(DBG,
"READ TPC data for sector %d (0x%x)",i+1,rdr);
96 LOG(DBG,
"Reading clusters");
98 sectorFirstHit[i+1] = nHits;
99 readITPCClustersFromEvpReader(rdr, i+1);
100 readClustersFromEvpReader(rdr, i+1);
104 LOG(DBG,
"Tracking...");
105 tracker->setClustersFromGl3Event(
this, i+1);
111 ret = tracker->processSector();
120 ret = tracker->fillTracks(szSECP_max, (
char *)sectp, 0);
125 int n = readSectorTracks((
char *)sectp);
129 LOG(WARN,
"Error reading tracker: sector %d\n",i,0,0,0,0);
137 finalizeReconstruction();
143 emc.readFromEvpReader(rdr, mem);
153 void gl3Event::readClustersFromEvpReader(
daqReader *rdr,
int sector)
160 dd = rdr->det(
"tpx")->get(
"legacy",sector);
162 dd = rdr->det(
"tpc")->get(
"legacy", sector);
167 pTPC = (
tpc_t *)dd->Void;
170 LOG(DBG,
"No data for sector %d check for TPC",sector);
174 LOG(DBG,
"have clusters? %p %d",pTPC, pTPC->has_clusters);
175 if(!pTPC->has_clusters)
return;
177 for(
int r=0;r<45;r++) {
179 for(
int j=0;j<pTPC->cl_counts[r];j++) {
180 tpc_cl *c = &pTPC->cl[r][j];
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;
197 gl3c->set(coordinateTransformer, sector, &sl3c);
209 int gl3Event::readITPCClustersFromEvpReader(
daqReader *rdr,
int sector) {
213 dd = rdr->det(
"itpc")->get(
"cld", sector);
215 while(dd->iterate()) {
218 int padrow = dd->row;
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;
231 gl3c->setITPCHit(coordinateTransformer, sec, padrow, pad, tb, charge, flags);
247 int mxHits,
int mxTracks)
248 : emc(inBemcCalib, inEemcCalib)
257 maxSectorNForTrackMerging = 1000000;
258 coordinateTransformer = inTrans;
260 setup( mxHits, mxTracks );
265 gl3Event::~gl3Event( )
267 if (
hit != 0 )
delete[]
hit ;
269 if ( trackContainer != 0 )
delete[] trackContainer;
270 if ( trackIndex != 0 )
delete[] trackIndex ;
279 gl3Track* gl3Event::getTrack (
int n ) {
280 if ( n < 0 || n > nTracks ) {
281 fprintf ( stderr,
" %d track index out of range \n", n );
288 gl3Hit* gl3Event::getHit (
int n ) {
289 if ( n < 0 || n > nHits ) {
290 fprintf ( stderr,
" %d hit index out of range \n", n );
296 gl3Sector* gl3Event::getSector (
int n ) {
297 if ( n < 0 || n > nSectors ) {
298 fprintf ( stderr,
" %d sector index out of range \n", n );
301 return &(sectorInfo[n]);
305 int gl3Event::getTrgCmd()
311 int gl3Event::getTrgWord()
313 return trgData.triggerWord;
316 int gl3Event::getZDC(
int n)
318 return trgData.ZDC[n];
321 int gl3Event::getCTB(
int n)
323 return trgData.CTB[n];
326 double gl3Event::getZDCVertex()
328 return ((
double)(trgData.ZDC[9] - trgData.ZDC[8]) + 21.3) * 3.3;
335 unsigned int gl3Event::getBXingLo()
337 return trgData.bunchXing_lo;
340 unsigned int gl3Event::getBXingHi()
342 return trgData.bunchXing_hi;
345 unsigned long long gl3Event::getBXing()
347 unsigned long long bx_hi_long = trgData.bunchXing_hi;
348 unsigned long long bx_lo_long = trgData.bunchXing_lo;
350 return (bx_hi_long << 32) | bx_lo_long;
363 int indexStore = -1 ;
366 for (
int i = 0 ; i < nTrk ; i++ ) {
367 lTrack->
set ( sector, trk ) ;
368 lTrack->id = sector * 1000 + abs(trk->id) ;
369 lTrack->para = ¶ ;
370 lTrack->sector = sector ;
376 if ( hitProcessing ) {
378 if ( abs(idTrack) < maxTracksSector )
379 indexStore = (sector-1)*maxTracksSector + abs(idTrack) ;
381 LOG(ERR,
" gl3Event::addTracks: max number of tracks per Sector reached %d reached", idTrack ,0,0,0,0) ;
388 if ( maxSectorNForTrackMerging > nTrk && idTrack < 0 ) {
390 fatterTrack = lTrack->merge ( trackContainer ) ;
392 if ( hitProcessing && indexStore > 0 ) {
393 trackIndex[indexStore] =
404 if ( hitProcessing && indexStore > 0 )
405 trackIndex[indexStore] = nTracks + 1;
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) ;
425 int gl3Event::fillTracks (
int maxBytes,
char* buffer,
unsigned int token ){
430 if ( nBytesNeeded > maxBytes ) {
431 LOG(ERR,
" gl3Event::writeTracks: %d bytes needed less than max = %d \n",
432 nBytesNeeded, maxBytes ,0,0,0) ;
439 head->xVert =
vertex.Getx();
440 head->yVert =
vertex.Gety();
441 head->zVert =
vertex.Getz();
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;
452 head->bh.length = (
sizeof(
struct L3_GTD)
461 for (
int i = 0 ; i < nTracks ; i++ ) {
462 if ( fabs(
track[i].z0) > 205 ) {
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 ;
491 head->nTracks = counter ;
495 return ((
char *)oTrack-buffer) ;
503 int gl3Event::readL3Data(
L3_P* header )
506 char* buffer = (
char *)header;
516 for ( i = 0 ; i < nSectors ; i++ ) {
517 length = header->sector[i].len ;
519 if ( length==0 ) continue ;
521 offset = 4 * header->sector[i].off ;
522 sectorP = (
L3_SECP *)&(buffer[offset]);
524 trackPointer = (
char *)sectorP + sectorP->trackp.off * 4 ;
527 int nSectorTracks = 0;
528 if (sectorP->trackp.off) {
529 nSectorTracks = readSectorTracks ( trackPointer ) ;
531 if ( nSectorTracks < 0 ) {
532 LOG(ERR,
"gl3Event:readEvent: error reading tracks, sector %d", i+1,0,0,0,0);
537 if ( hitProcessing && sectorP->sl3clusterp.off ) {
538 hitPointer = (
char *)sectorP + sectorP->sl3clusterp.off * 4 ;
539 readSectorHits ( hitPointer, nSectorTracks ) ;
544 if(header->bh.format_number>=5 && header->trig.len){
547 trgData.readL3P(header);
552 if( header->bh.format_number>=7 ){
554 emc.readRawData(header);
567 for (
int i = 0 ; i < nTracks ; i++ ) {
569 radius = coordinateTransformer->
570 GetRadialDistanceAtRow(
track[i].innerMostRow-1) ;
572 track[i].updateToRadius ( radius ) ;
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",
595 int gl3Event::finalizeReconstruction()
599 if (vertexFinder & 0x01)
604 if ((vertexFinder & 0x02) && lmv) {
606 lmv->makeVertex(
this);
609 lmVertex.Setx(vtx.x);
610 lmVertex.Sety(vtx.y);
611 lmVertex.Setz(vtx.z);
616 vertex_ftf.x =
vertex.Getx();
617 vertex_ftf.y =
vertex.Gety();
618 vertex_ftf.z =
vertex.Getz();
621 for (
int i=0 ; i<getNTracks() ; i++) {
622 getTrack(i)->setDca(vertex_ftf);
634 int gl3Event::readSectorHits (
char* buffer,
int nSectorTracks ){
639 if ( !coordinateTransformer ) {
640 LOG(ERR,
"gl3Event::readSectorHits: there is not Coordinate Transformer",0,0,0,0,0);
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 ) ;
654 int sector = head->bh.bank_id;
655 int nSectorHits = head->nrClusters_in_sector ;
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) ;
671 for (
int i = 0 ; i < nSectorHits ; i++ ) {
672 hitP = &(cluster[i]) ;
676 if ( hitProcessing > 1 ) {
677 gHitP = &(
hit[nHits+i]);
678 gHitP->set (coordinateTransformer, sector, hitP);
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) ;
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) ;
696 int index = trackIndex[indexStore] - 1 ;
697 if ( index < 0 || index > nTracks ) continue ;
702 if ( hitProcessing > 1 ) {
703 if (
track[index].firstHit == 0 )
704 track[index].firstHit = (
void *)gHitP ;
706 ((
gl3Hit *)(
track[index].lastHit))->nextHit = (
void *)gHitP ;
707 track[index].lastHit = (
void *)gHitP ;
708 gHitP->trackId =
track[index].id ;
713 hitP->trackId =
track[index].id ;
717 nHits += nSectorHits ;
727 int gl3Event::readSectorTracks (
char* buffer ){
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 ) ;
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) ;
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 ;
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 ;
767 char* pointer = head->banks[0].off * 4 + buffer ;
770 if ( (head->banks[0].len > 0) && (head->bh.format_number > 0) ) {
772 nSectorTracks = (4 * head->banks[0].len -
sizeof(
struct bankHeader))
775 else nSectorTracks = 0 ;
779 if ( nSectorTracks > 0 ) {
782 addTracks ( sector, nSectorTracks, localTrack ) ;
786 return sectorP->nTracks ;
793 int gl3Event::makeVertex (){
807 vertex.Setxyz(0.0,0.0,0.0);
810 for(
int iter = 0 ; iter<2; iter++ ) {
812 for(
int trkcnt = 0 ; trkcnt<getNTracks(); trkcnt++ ) {
813 gTrack = getTrack(trkcnt);
818 if ( gTrack->nHits > minNoOfHitsOnTrackUsedForVertexCalc &&
819 gTrack->pt > minMomUsedForVertexCalc) {
821 closestHit = gTrack->closestApproach(getVertex().Getx(),
827 hVz->
Fill(closestHit.z,1.0);
828 hVx->
Fill(closestHit.x,1.0);
829 hVy->
Fill(closestHit.y,1.0);
834 vertex.Setxyz(hVx->getWeightedMean(6.0),
835 hVy->getWeightedMean(6.0),
836 hVz->getWeightedMean(4.0));
847 int gl3Event::resetEvent ( ){
851 nMergableTracks = 0 ;
855 memset(sectorFirstHit, 0,
sizeof(sectorFirstHit));
858 memset(trackContainer, 0,
859 para.nPhiTrackPlusOne*para.nEtaTrackPlusOne*
sizeof(
FtfContainer));
862 if ( hitProcessing ) {
863 memset ( trackIndex, 0, maxTracksSector*nSectors*
sizeof(
int) ) ;
869 for (
int i=0; i<16; i++)
872 for (
int i=0; i<240; i++)
884 vertex.Setxyz(0.0, 0.0, 0.0);
894 int gl3Event::setup (
int mxHits,
int mxTracks )
897 if ( mxHits < 0 || mxHits > 1000000 ) {
898 LOG(ERR,
" gl3Event::setup: maxHits %d out of range \n", maxHits,0,0,0,0 ) ;
902 if ( mxTracks < 0 || mxTracks > 1000000 ) {
903 LOG(ERR,
" gl3Event::setup: maxTracks %d out of range \n", maxTracks,0,0,0,0 );
909 maxTracks = mxTracks ;
910 maxTracksSector = maxTracks*2/ nSectors ;
913 trackIndex =
new int[maxTracksSector*nSectors];
919 para.nPhiTrackPlusOne = para.nPhiTrack + 1 ;
920 para.nEtaTrackPlusOne = para.nEtaTrack + 1 ;
924 para.phiSliceTrack = (para.phiMaxTrack - para.phiMinTrack)/para.nPhiTrack;
925 para.etaSliceTrack = (para.etaMaxTrack - para.etaMinTrack)/para.nEtaTrack;
927 int nTrackVolumes = para.nPhiTrackPlusOne* para.nEtaTrackPlusOne ;
929 if(trackContainer == NULL) {
930 LOG(ERR,
"Problem with memory allocation... exiting\n",0,0,0,0,0) ;
934 para.ptMinHelixFit = 1.e60 ;
941 minNoOfHitsOnTrackUsedForVertexCalc=14;
942 minMomUsedForVertexCalc=0.25;
947 strcpy ( hid,
"Vertex_Vz" ) ;
948 strcpy ( title,
"Vertex_Vz" ) ;
949 hVz =
new gl3Histo ( hid, title, 400, -200., 200. ) ;
951 strcpy ( hid,
"Vertex_Vx" ) ;
952 strcpy ( title,
"Vertex_Vx" ) ;
953 hVx =
new gl3Histo ( hid, title, 100,-10,10);
955 strcpy ( hid,
"Vertex_Vy" ) ;
956 strcpy ( title,
"Vertex_Vy" ) ;
957 hVy =
new gl3Histo ( hid, title, 100,-10,10);
void set(short sectorIn, local_track *trk)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
int Reset()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
int readFromEvpReader(daqReader *rdr, float bField=1000)
void addTracks(bool cuts=false)
Add tracks to the list of the rendered objects from current MuDst event.
int Fill(double x, double weight=1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ...