GltValsTool class Reference

00253 {
00254     StatusCode sc = StatusCode::SUCCESS;
00255 
00256     // m_pEventSvc already checked by doCalcIfNotDone, no need to repeat
00257 
00258     int nLayers  = m_tkrGeom->numLayers();
00259     int nXTowers = m_tkrGeom->numXTowers();
00260     int nYTowers = m_tkrGeom->numYTowers();
00261     int nTowers  = nXTowers*nYTowers;
00262 
00263 
00264     int iTrig_tower  = -1;
00265     int iTrig_layer  = -1;
00266     int iTrig_xTower = -1;
00267     int iTrig_yTower = -1; 
00268     int iTrig_type   =  0;
00269 
00270     Trig_layer  = iTrig_layer;
00271     Trig_tower  = iTrig_tower;
00272     Trig_xTower = iTrig_xTower;
00273     Trig_yTower = iTrig_yTower;
00274     Trig_type   = iTrig_type;
00275 
00276     // Recover EventHeader Pointer
00277     SmartDataPtr<Event::EventHeader> 
00278         pEvent(m_pEventSvc, EventModel::EventHeader);
00279     unsigned int word  = 0;
00280     unsigned int word2 = _unset;
00281     
00282     if(pEvent) {
00283     word  = pEvent->trigger();
00284     word2 = pEvent->triggerWordTwo();
00285      // by default, engines are initialized to "_unset"; see zeroVals() above
00286     Trig_gemprescale     = pEvent->gemPrescale();
00287     Trig_gltprescale     = pEvent->gltPrescale();
00288     Trig_prescaleexpired = pEvent->prescaleExpired();
00289     }
00290 
00291     // This is the same as the old GltWord
00292     // actually only 6 bits, but no harm (I think!)  
00293     // note that the trigger word from the header has 3 GEM_offset fields that we unpack here:
00294     Trig_word = word & enums::GEM_mask;  // the GltWord, set for simulation from hits or digis
00295 
00296     Trig_GemSummary = (word >> enums::GEM_offset) & enums::GEM_mask; // the GEM condition word, same as previous if simulation
00297     
00298     Trig_engine = (word >> (2*enums::GEM_offset)) & enums::GEM_mask; // the trigger engine number
00299 
00300     unsigned int Trig_gltengine = word2 & enums::ENGINE_mask;
00301     // If the engine number is set, then we will store this one from the TrgConfigSvc in the tuple
00302     if (Trig_gltengine != enums::ENGINE_unset) Trig_engine = Trig_gltengine;
00303     Trig_gemengine = ((word2 >> enums::ENGINE_offset) & enums::ENGINE_mask); // GEM trigger engine number
00304 
00305     SmartDataPtr<LsfEvent::MetaEvent> 
00306         metaTds(m_pEventSvc, "/Event/MetaEvent");
00307     if(metaTds) {
00308         Trig_sourcegps =  metaTds->time().current().sourceGps();
00309         Trig_compressedEventSize = metaTds->compressedSize();
00310     }
00311 
00312     SmartDataPtr<LdfEvent::EventSummaryData> 
00313         eventSummary(m_pEventSvc, "/Event/EventSummary"); 
00314     if(eventSummary) {
00315         Trig_evtFlags =  eventSummary->eventFlags();
00316         Trig_eventSize = eventSummary->eventSizeInBytes();
00317     }
00318 
00319     // GemDeltaTimes
00320     SmartDataPtr<LdfEvent::Gem> gemTds(m_pEventSvc, "/Event/Gem");
00321     if(gemTds) {
00322         Trig_gemDeltaEventTime      = gemTds->deltaEventTime();
00323         Trig_gemDeltaWindowOpenTime = gemTds->deltaWindowOpenTime();
00324     }
00325 
00326     SmartDataPtr<Event::TkrClusterCol>   
00327         pClusters(m_pEventSvc,EventModel::TkrRecon::TkrClusterCol);
00328 
00329     // everything from here on out uses Clusters
00330     if(!pClusters) return sc;
00331 
00332     bool three_in_a_row = ((word & enums::b_Track)!=0);
00333 
00334     int tower, layer;
00335     // needs to be recast into an indexed vector maybe
00336     short x_hits[_nTowers][_nLayers]; 
00337     short y_hits[_nTowers][_nLayers]; 
00338     for(tower=0; tower<_nTowers; ++tower) {
00339         for(layer=0; layer<_nLayers; ++layer){
00340             x_hits[tower][layer]=0;
00341             y_hits[tower][layer]=0;
00342         }
00343     }
00344 
00345     //Search for x-y paired layer.... 
00346     if (pClusters && three_in_a_row)
00347     {
00348         // Make a hit count per plane, per tower.  
00349         layer = nLayers;
00350         while(layer--)
00351         {
00352             Event::TkrClusterVec xHitList = 
00353                 m_clusTool->getClusters(idents::TkrId::eMeasureX,layer);
00354             Event::TkrClusterVec yHitList = 
00355                 m_clusTool->getClusters(idents::TkrId::eMeasureY,layer);
00356 
00357             int x_hitCount = xHitList.size(); 
00358             int y_hitCount = yHitList.size();
00359             if(x_hitCount > 0 && y_hitCount > 0) {
00360 
00361                 int hit;
00362                 // x hits
00363                 for(hit=0; hit<x_hitCount; ++hit){
00364                     tower = xHitList[hit]->tower();
00365                     x_hits[tower][layer]++; 
00366                 }
00367                 // y hits
00368                 for(hit=0; hit<y_hitCount; ++hit){
00369                     tower = yHitList[hit]->tower();
00370                     y_hits[tower][layer]++;
00371                 }
00372             }
00373         }
00374 
00375         // Now search for the 3-in-a-rows (need to look for x & y hits in each...)
00376         // If more than one tower triggers, use the one with the highest layer number
00377         int nTotTrigs = 0;   // total number of possible triggers
00378         int nTowerTrigs = 0; // total number of towers that could trigger
00379         for(tower = 0; tower<nTowers; ++tower) {
00380             int nTrigsThis = 0;
00381             for(layer=0; layer<nLayers-2; ++layer){
00382                 if( (x_hits[tower][layer]   > 0 && y_hits[tower][layer]   > 0 ) &&
00383                     (x_hits[tower][layer+1] > 0 && y_hits[tower][layer+1] > 0 ) &&
00384                     (x_hits[tower][layer+2] > 0 && y_hits[tower][layer+2] > 0 )) 
00385                 {
00386                     nTrigsThis++;
00387                     if(layer > iTrig_layer) {
00388                         iTrig_layer = layer;
00389                         iTrig_tower = tower;
00390                     }
00391                 }
00392             }
00393             nTotTrigs += nTrigsThis;
00394             if(nTrigsThis>0) nTowerTrigs++;
00395         }
00396         Trig_total  = nTotTrigs;
00397         Trig_numTowers = nTowerTrigs;
00398 
00399         // Now classify according to tower type
00400         // new classification is number of exposed sides
00401         if(iTrig_tower >= 0) {
00402             iTrig_type = m_tkrGeom->getTowerType(iTrig_tower);
00403             idents::TowerId towerId = idents::TowerId(iTrig_tower);
00404             iTrig_xTower = towerId.ix();
00405             iTrig_yTower = towerId.iy();
00406         }
00407 
00408         // Now find the average location of all hits
00409         // This is probably too coarse to be useful
00410         double x_sum, y_sum, z_sum, wts, g1, g2, L11, L12, L22, zmax, zmin;
00411         x_sum = y_sum = z_sum = wts = g1 = g2 = L11 = L12 = L22 = 0.;
00412         zmax = 0.;
00413         zmin = 20.; 
00414         CLHEP::HepSymMatrix moments(3,0); 
00415 
00416         if(iTrig_tower >= 0) {
00417             for(tower = 0; tower<_nTowers; ++tower) {
00418                 idents::TowerId towerId = idents::TowerId(tower);
00419                 int ix = towerId.ix();
00420                 int iy = towerId.iy();
00421                 if (ix+1==nXTowers || iy+1==nYTowers) continue;
00422                 for(layer=0; layer<nLayers; ++layer){
00423                     int x_counts = x_hits[tower][layer];
00424                     int y_counts = y_hits[tower][layer];
00425                     double weight = x_counts+y_counts;
00426                     if(x_counts>0 && y_counts>0) { // Don't count random noise hits
00427                         wts += weight;
00428                         double x = (ix+1)*5.;
00429                         double y = (iy+1)*5.;
00430                         double z = (layer+1);
00431                         x_sum += x * weight;
00432                         y_sum += y * weight;
00433                         z_sum += z * weight;
00434                         g1    += 1.;
00435                         g2    += z;
00436                         L11   += 1./weight;
00437                         L12   += z /weight;
00438                         L22   += z*z/weight; 
00439                         if(zmax < z ) zmax = z;
00440                         if(zmin > z ) zmin = z; 
00441 
00442                     }
00443                 }
00444             }
00445 
00446             // Now form moment tensor, etc..... 
00447             if(wts > 0) {
00448                 x_sum /= wts;
00449                 y_sum /= wts;
00450                 z_sum /= wts;
00451                 for(tower = 0; tower<_nTowers; ++tower) {
00452                     idents::TowerId towerId = idents::TowerId(tower);
00453                     int ix = towerId.ix();
00454                     int iy = towerId.iy();
00455                     if (ix+1==nXTowers || iy+1==nYTowers) continue;
00456                     for(layer=0; layer<nLayers-2; ++layer){
00457                         int x_counts = x_hits[tower][layer];
00458                         int y_counts = y_hits[tower][layer];
00459                         double weight = x_counts+y_counts;
00460                         if(x_counts>0 && y_counts>0) { // Don't count random noise hits
00461                             double xres = (ix+1)*5.  - x_sum;
00462                             double yres = (iy+1)*5.  - y_sum;
00463                             double zres = (layer+1)  - z_sum;
00464                             moments(1,1) += xres * xres * weight;
00465                             moments(1,2) += xres * yres * weight;
00466                             moments(1,3) += xres * zres * weight;
00467                             moments(2,2) += yres * yres * weight;
00468                             moments(2,3) += yres * zres * weight;
00469                             moments(3,3) += zres * zres * weight;
00470                         }
00471                     }
00472                 }
00473                 moments /= wts;
00474 
00475                 // Now find transformation to diagonalize
00476                 CLHEP::HepMatrix U(3,3);
00477                 U = diagonalize(&moments);
00478 
00479                 // Now get directions, etc. 
00480                 double Det = L11*L22 - L12*L12;
00481                 if(fabs(Det) > 0) {
00482                     double z_hit_slope = (g2*L11 - g1*L12)/Det;
00483                     Trig_zDir = z_hit_slope*(zmax-zmin)*wts; 
00484                 }
00485                 int sm = 1;
00486                 if(moments(2,2) > moments(1,1))   sm = 2;
00487                 if(moments(3,3) > moments(sm,sm)) sm = 3;
00488 
00489                 Trig_moment = moments(sm,sm);
00490                 //            Vector t_moment(U(1,sm),U(2,sm),U(3,sm));
00491             }
00492         }
00493     }
00494 
00495     Trig_layer  = iTrig_layer;
00496     Trig_tower  = iTrig_tower;
00497     Trig_xTower = iTrig_xTower;
00498     Trig_yTower = iTrig_yTower;
00499     Trig_type   = iTrig_type;
00500 
00501     return sc;
00502 }
GltValsTool Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	GltValsTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual	~GltValsTool ()
StatusCode	initialize ()
StatusCode	calculate ()
	calculate all values; implemented by each XxxValsTool
void	zeroVals ()
	clear map values
Private Attributes
ITkrGeometrySvc *	m_tkrGeom
	pointer to tracker geometry
float	Trig_word
float	Trig_GemSummary
float	Trig_evtFlags
float	Trig_tower
float	Trig_xTower
float	Trig_yTower
float	Trig_layer
float	Trig_total
float	Trig_numTowers
float	Trig_type
float	Trig_moment
float	Trig_zDir
int	Trig_engine
int	Trig_gemengine
int	Trig_gltprescale
int	Trig_gemprescale
int	Trig_prescaleexpired
int	Trig_sourcegps
int	Trig_gemDeltaEventTime
int	Trig_gemDeltaWindowOpenTime
int	Trig_eventSize
int	Trig_compressedEventSize
ITkrQueryClustersTool *	m_clusTool