10 #include "StEEmcPointMap.h"
14 #include "StRoot/St_base/StMessMgr.h"
16 #include "StEEmcStripEndPointData.h"
17 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
18 #include "StRoot/StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
19 #include "StRoot/StEEmcPool/StEEmcGeoId/StEEmcGeoId.h"
22 #ifndef PI_OVER_TWELVE
23 #define PI_OVER_TWELVE 0.2617993877991494263
27 #ifndef TWO_PI_OVER_TWELVE
28 #define TWO_PI_OVER_TWELVE 0.5235987755982988527
33 #define EDGE_FACTOR 1.025
37 void StEEmcPointMap_t::getStripLineParam( Float_t stripGeoId, Float_t& a, Float_t& b ) {
38 Float_t x1, x2, y1, y2;
39 getStripEndPoints( stripGeoId, x1, y1, x2, y2 );
41 Float_t theta = atan2( y2 - y1, x2 - x1 );
45 theta = 0.0872664*
static_cast< Int_t
>( theta/0.0872664 + 0.5 );
52 void StEEmcPointMap_t::getStripEndPoints( Float_t stripGeoId, Float_t& x1, Float_t& y1, Float_t& x2, Float_t& y2 ) {
55 if( mStripDataVec.empty() )
59 Short_t stripSector, stripIndex;
63 StEEmcGeoId_t::decodeSmd( stripGeoId, stripSector, layerIsV, stripIndex );
65 if( (UInt_t)stripIndex > mStripDataVec.size() ){
66 LOG_ERROR <<
"Asked for invalid strip: " << endm;
67 LOG_ERROR <<
"\tGeoId " << stripGeoId <<
", sector " << stripSector <<
" layer " << ( layerIsV ?
'V' :
'U' ) <<
", index " << stripIndex << endm;
68 LOG_ERROR <<
"\tStrip end point data currently loaded for " << mStripDataVec.size() <<
" strips." << endm;
69 x1 = y1 = x2 = y2 = 0;
74 StEEmcStripEndPointDataVec_t::const_iterator stripDataSetIter = mStripDataVec.begin();
75 std::advance( stripDataSetIter, stripIndex );
79 Float_t dphiU = -(stripSector * TWO_PI_OVER_TWELVE - PI_OVER_TWELVE);
81 Float_t x1_ = stripDataSetIter->getX();
82 Float_t y1_ = stripDataSetIter->getY1();
83 Float_t x2_ = stripDataSetIter->getX();
84 Float_t y2_ = stripDataSetIter->getY2();
86 Float_t frac = stripGeoId - (Int_t)(stripGeoId);
89 if( stripDataSetIter != mStripDataVec.end() ){
90 Float_t x1b_ = stripDataSetIter->getX();
91 Float_t y1b_ = stripDataSetIter->getY1();
92 Float_t x2b_ = stripDataSetIter->getX();
93 Float_t y2b_ = stripDataSetIter->getY2();
95 x1_ += frac*(x1b_ - x1_);
96 y1_ += frac*(y1b_ - y1_);
97 x2_ += frac*(x2b_ - x2_);
98 y2_ += frac*(y2b_ - y2_);
112 x1 = x1_*cos( dphiU ) - y1_*sin( dphiU );
113 y1 = y1_*cos( dphiU ) + x1_*sin( dphiU );
115 x2 = x2_*cos( dphiU ) - y2_*sin( dphiU );
116 y2 = y2_*cos( dphiU ) + x2_*sin( dphiU );
120 Float_t StEEmcPointMap_t::getStripDCA( Float_t pointX, Float_t pointY, Short_t stripGeoId )
const {
126 Short_t stripSector, stripIndex;
130 StEEmcGeoId_t::decodeSmd( stripGeoId, stripSector, layerIsV, stripIndex );
133 StEEmcStripEndPointDataVec_t::const_iterator stripDataSetIter = mStripDataVec.begin();
134 std::advance( stripDataSetIter, stripIndex );
141 Float_t dphiU = stripSector * TWO_PI_OVER_TWELVE - PI_OVER_TWELVE;
142 Float_t pointXtrans = pointX*cos( dphiU ) - pointY*sin( dphiU );
143 Float_t pointYtrans = pointY*cos( dphiU ) + pointX*sin( dphiU );
146 Float_t temp = pointXtrans;
147 pointXtrans = pointYtrans;
157 if( pointYtrans > stripDataSetIter->getY2() ){
159 dX = pointXtrans - stripDataSetIter->getX();
160 dY = pointYtrans - stripDataSetIter->getY2();
162 }
else if( pointYtrans < stripDataSetIter->getY1() ){
164 dX = pointXtrans - stripDataSetIter->getX();
165 dY = pointYtrans - stripDataSetIter->getY1();
170 dX = pointXtrans - stripDataSetIter->getX();
175 Bool_t isEdge = ( !layerIsV && (stripSector == 2 || stripSector == 8 )) || (layerIsV && (stripSector == 3 || stripSector == 9 ));
176 Float_t phi = atan2( pointYtrans, pointXtrans );
179 if( isEdge && phi < TWO_PI_OVER_TWELVE*EDGE_FACTOR ){
180 Float_t altY1 = stripDataSetIter->getX()*tan( TWO_PI_OVER_TWELVE*EDGE_FACTOR );
185 if( altY1 > stripDataSetIter->getY1() )
186 dY = pointYtrans-altY1;
189 return sqrt( dX*dX + dY*dY );
193 Float_t StEEmcPointMap_t::getTowerDCA( Float_t x, Float_t y, Short_t towerGeoId )
const {
198 StEEmcGeoId_t::decodeTow( towerGeoId, sector, subsector, etaBin );
202 Float_t dX = x - towerCenter.X();
203 Float_t dY = y - towerCenter.Y();
205 return sqrt( dX*dX + dY*dY );
209 Short_t StEEmcPointMap_t::getStripsNearestPoint( Float_t x, Float_t y, Float_t r )
const{
212 std::set< Short_t > stripSet;
214 for( ; dist < r && geoId < 0; dist += 0.5 ){
215 getStripsNearPoint( x, y, dist, stripSet );
217 if( stripSet.size() == 1 ){
218 geoId = *(stripSet.begin());
219 }
else if (!stripSet.empty() ){
220 Float_t smallestDist = 2*r;
221 Short_t geoIdOfSmallest = -1;
223 std::set< Short_t >::iterator iter;
224 for( iter = stripSet.begin(); iter != stripSet.end(); ++iter ){
225 Float_t thisDist = getStripDCA( x, y, *iter );
226 if( thisDist < smallestDist ){
227 smallestDist = thisDist;
228 geoIdOfSmallest = *iter;
231 geoId = geoIdOfSmallest;
239 void StEEmcPointMap_t::getStripsNearPoint( Float_t x, Float_t y, Float_t r, std::set< Short_t >& stripSet )
const {
241 Float_t phi = atan2( y, x );
244 Float_t sectorFrac = ( 2.5 - phi/TWO_PI_OVER_TWELVE );
245 Short_t sector = (Short_t)sectorFrac;
248 if( sectorFrac <= 0 )
251 sector %= kEEmcNumSectors;
253 sector += kEEmcNumSectors;
256 Float_t dphiU = sector * TWO_PI_OVER_TWELVE - PI_OVER_TWELVE;
257 Float_t x2 = x*cos( dphiU ) - y*sin( dphiU );
258 Float_t y2 = y*cos( dphiU ) + x*sin( dphiU );
261 addSmdContribution( x2, y2, r, sector,
'U', 0, stripSet );
262 addSmdContribution( x2, y2, r, sector,
'V', 0, stripSet );
264 if( mCrossSectorBoundaries ){
269 Float_t x3 = x2*cos( TWO_PI_OVER_TWELVE ) - y2*sin( TWO_PI_OVER_TWELVE );
270 Float_t y3 = y2*cos( TWO_PI_OVER_TWELVE ) + x2*sin( TWO_PI_OVER_TWELVE );
273 Short_t newSector = (sector + 1) % kEEmcNumSectors;
274 addSmdContribution( x3, y3, r, newSector,
'U', 1, stripSet );
275 addSmdContribution( x3, y3, r, newSector,
'V', 1, stripSet );
278 Float_t x4 = x2*cos( -TWO_PI_OVER_TWELVE ) - y2*sin( -TWO_PI_OVER_TWELVE );
279 Float_t y4 = y2*cos( -TWO_PI_OVER_TWELVE ) + x2*sin( -TWO_PI_OVER_TWELVE );
283 newSector = (sector - 1 + kEEmcNumSectors) % kEEmcNumSectors;
284 addSmdContribution( x4, y4, r, newSector,
'U', -1, stripSet );
285 addSmdContribution( x4, y4, r, newSector,
'V', -1, stripSet );
295 void StEEmcPointMap_t::addSmdContribution( Float_t x, Float_t y, Float_t r, Short_t sector, Char_t layer, Short_t sectorSide, std::set< Short_t >& stripSet )
const {
296 std::set< Short_t > stripIdxSet;
297 std::set< Short_t >::iterator stripIdxSetIter;
307 StEEmcStripEndPointDataVec_t::const_iterator low_end = std::lower_bound( mStripDataVec.begin(), mStripDataVec.end(),
StEEmcStripEndPointData_t( 0, x-r, 0, 0 ) );
308 StEEmcStripEndPointDataVec_t::const_iterator up_end = std::upper_bound( mStripDataVec.begin(), mStripDataVec.end(),
StEEmcStripEndPointData_t( 0, x+r, 0, 0 ) );
309 StEEmcStripEndPointDataVec_t::const_iterator stripDataSetIter;
311 Bool_t isEdge = (layer ==
'U' && (sector == 2 || sector == 8 )) || (layer ==
'V' && (sector == 3 || sector == 9 ));
313 for( ; low_end != up_end; ++low_end )
314 stripIdxSet.insert( low_end->getStripIndex() );
317 stripIdxSet.erase( kEEmcNumStrips-1 );
321 stripIdxSetIter = stripIdxSet.lower_bound( kEEmcNumEdgeStrips-3 );
323 if( stripIdxSetIter != stripIdxSet.end() )
324 stripIdxSet.erase( stripIdxSetIter, stripIdxSet.end() );
328 Float_t thresRsquared = r*r;
330 Float_t phi = atan2( y, x );
333 for( stripIdxSetIter = stripIdxSet.begin(); stripIdxSetIter != stripIdxSet.end(); ++stripIdxSetIter ){
339 stripDataSetIter = mStripDataVec.begin();
340 std::advance( stripDataSetIter, *stripIdxSetIter );
343 Float_t xdist = x - stripDataSetIter->getX();
345 if( y < stripDataSetIter->getY1() )
346 ydist = y - stripDataSetIter->getY1();
347 if( y > stripDataSetIter->getY2() )
348 ydist = y - stripDataSetIter->getY2();
352 if( isEdge && phi < TWO_PI_OVER_TWELVE*EDGE_FACTOR ){
353 Float_t altY1 = stripDataSetIter->getX()*tan( TWO_PI_OVER_TWELVE*EDGE_FACTOR );
361 if( altY1 > stripDataSetIter->getY1() )
365 Float_t dist_sq = xdist*xdist + ydist*ydist;
367 if( dist_sq < thresRsquared ){
368 stripSet.insert( StEEmcGeoId_t::encodeSmd( sector, (layer ==
'V'), *stripIdxSetIter ) );
370 stripDataSetIter = mStripDataVec.begin();
371 std::advance( stripDataSetIter, *stripIdxSetIter );
377 void StEEmcPointMap_t::getTowersNearPoint( Float_t x, Float_t y, Float_t r, std::set< Short_t >& towerSet ){
380 Float_t phi = atan2( y, x );
382 Float_t rOfXY = sqrt( x*x + y*y );
385 LOG_ERROR <<
"Invalid point, distance combination: (" << x <<
", " << y <<
"), " << r << endm;
390 Float_t dphi = asin( r/rOfXY );
392 Short_t kNumPhiBins = kEEmcNumSectors*kEEmcNumSubSectors;
395 Float_t phiFrac1 = ( 12.5 - (phi+dphi)/TWO_PI_OVER_TWELVE*kEEmcNumSubSectors );
396 Short_t phiBin1 = (Short_t)phiFrac1;
403 phiBin1 += kNumPhiBins;
404 phiBin1 %= kNumPhiBins;
406 Float_t phiFrac2 = ( 12.5 - (phi-dphi)/TWO_PI_OVER_TWELVE*kEEmcNumSubSectors );
407 Short_t phiBin2 = (Short_t)phiFrac2;
414 phiBin2 += kNumPhiBins;
415 phiBin2 %= kNumPhiBins;
417 if( phiBin2 < phiBin1 )
418 phiBin2 += kNumPhiBins;
421 Float_t minTheta = atan2( rOfXY - r, kEEmcZSMD );
422 Float_t maxTheta = atan2( rOfXY + r, kEEmcZSMD );
424 Float_t minEta = - log( tan( maxTheta / 2 ));
425 Float_t maxEta = - log( tan( minTheta / 2 ));
428 Short_t minEtaBin = 0;
429 Short_t maxEtaBin = kEEmcNumEtas-1;
431 if( minEta > etaBinRange[ kEEmcNumEtas ] && maxEta < etaBinRange[ 0 ] ){
435 for( ; minEtaBin <= kEEmcNumEtas+1 && etaBinRange[ minEtaBin ] > maxEta; ++minEtaBin );
441 for( maxEtaBin = minEtaBin; maxEtaBin <= kEEmcNumEtas+1 && etaBinRange[ maxEtaBin ] > minEta; ++maxEtaBin );
449 for( Short_t etaBin = minEtaBin; etaBin <= maxEtaBin; ++etaBin )
450 for( Short_t phiBin = phiBin1; phiBin <= phiBin2; ++phiBin )
451 towerSet.insert( StEEmcGeoId_t::encodeTow( phiBin%kNumPhiBins, etaBin ) );
457 Short_t StEEmcPointMap_t::getTowerNearestPoint( Float_t x, Float_t y ){
460 Float_t phi = atan2( y, x );
461 Float_t rOfXY = sqrt( x*x + y*y );
463 Short_t kNumPhiBins = kEEmcNumSectors*kEEmcNumSubSectors;
466 Float_t phiFrac = ( 12.5 - phi/TWO_PI_OVER_TWELVE*kEEmcNumSubSectors );
467 Short_t phiBin = (Short_t)phiFrac;
474 phiBin += kNumPhiBins;
475 phiBin %= kNumPhiBins;
478 Float_t theta = atan2( rOfXY, kEEmcZSMD );
479 Float_t eta = - log( tan( theta / 2 ));
484 if( eta > etaBinRange[ kEEmcNumEtas ] && eta < etaBinRange[ 0 ] ){
487 for( ; etaBin < kEEmcNumEtas && eta < etaBinRange[ etaBin+1 ]; ++etaBin );
489 towIdx = StEEmcGeoId_t::encodeTow( phiBin%kNumPhiBins, etaBin );
495 Short_t StEEmcPointMap_t::getSectorOfPoint( Float_t x, Float_t y ){
496 Float_t phi = atan2( y, x );
498 Float_t sectorFrac = ( 2.5 - phi/TWO_PI_OVER_TWELVE );
499 Short_t sector = (Short_t)sectorFrac;
502 if( sectorFrac <= 0 )
505 sector %= kEEmcNumSectors;
507 sector += kEEmcNumSectors;
512 void StEEmcPointMap_t::convertStripUVtoXY( Short_t sector, Float_t u, Float_t v, Float_t& x, Float_t& y ) {
516 Float_t uGeoId = StEEmcGeoId_t::encodeSmd( sector, 0, (Int_t)(u) ) + u - (Int_t)(u);
517 Float_t vGeoId = StEEmcGeoId_t::encodeSmd( sector, 1, (Int_t)(v) ) + v - (Int_t)(v);
520 Float_t uA, uB, vA, vB;
521 getStripLineParam( uGeoId, uA, uB );
522 getStripLineParam( vGeoId, vA, vB );
525 x = (uB - vB)/(vA - uA);
535 const StEEmcStripEndPointDataVec_t StEEmcPointMap_t::mStripDataVec;
538 StEEmcPointMap_t::StEEmcPointMap_t( Bool_t crossSectorBoundaries ) : mCrossSectorBoundaries( crossSectorBoundaries ) {
543 void StEEmcPointMap_t::loadData(){
544 if( mStripDataVec.empty() ){
546 StEEmcStripEndPointDataVec_t *stripDataSetPtr =
const_cast< StEEmcStripEndPointDataVec_t*
>( &mStripDataVec );
837 std::sort( stripDataSetPtr->begin(), stripDataSetPtr->end() );
TVector3 getTowerCenter(const UInt_t sec, const UInt_t sub, const UInt_t etabin) const
static EEmcGeomSimple & Instance()
returns a reference to a static instance of EEmcGeomSimple