StRoot: StEEmcUtil/StEEmcSmd/EEmcSmdGeom.cxx Source File

00001 
00035 #include "Stiostream.h"
00036 #include "EEmcSmdGeom.h"
00037 #include "EEmcStripGeom.h"
00038 #include <assert.h>
00039 #include <TMath.h>
00040 
00041 // decouple from StarClassLibrary
00042 //#include "PhysicalConstants.h"  
00043 #ifndef HEP_SYSTEM_OF_UNITS_H
00044 static const double     radian      = 1.;
00045 static const double     pi          = M_PI; // from <math.h>
00046 static const double     degree      = (M_PI/180.0)*radian;
00047 #endif
00048 
00049 
00050 #include "StMessMgr.h"
00051 
00052 
00054 ClassImp(EEmcSmdGeom)
00055 
00056 EEmcSmdGeom::EEmcSmdGeom()
00057     : TObject()
00058 { 
00059   for(int iSec = 0; iSec < kEEmcNumSectors; iSec++) mIsSectorIn[iSec] = true;
00060 };
00062 EEmcSmdGeom::~EEmcSmdGeom(){ 
00063   delete sInstance;
00064   sInstance = 0;
00065 }
00066 
00068 void EEmcSmdGeom::init(){ 
00069         buildSmdGeom(); 
00070 }
00071 
00072 EEmcSmdGeom* EEmcSmdGeom::sInstance = 0;        
00073 // all setctors
00074 EEmcSmdGeom* EEmcSmdGeom::instance() {  
00075     if(!sInstance){
00076         sInstance = new EEmcSmdGeom();
00077         sInstance->init();
00078    }
00079    return sInstance;
00080 }     
00081 
00082 // selected sectors
00083 EEmcSmdGeom* EEmcSmdGeom::instance(intVec sectorIdVec) {        
00084     if(!sInstance){
00085         sInstance = new EEmcSmdGeom();
00086         sInstance->setSectors(sectorIdVec);    
00087         sInstance->init();
00088    }
00089    return sInstance;
00090 }     
00091 
00092 // build a glabal geometry database from local coordinates
00093 void EEmcSmdGeom::buildSmdGeom(){
00094   mEEmcSmdParam.stripWidth = 0.5;
00095   for(int iPlane=0; iPlane<kEEmcNumSmdPlanes; iPlane++) 
00096                  mEEmcSmdParam.rOffset[iPlane] = kEEmcSmdROffset[iPlane];
00097 
00098   float x0[kEEmcNumStrips];
00099   float y1[kEEmcNumStrips];
00100   float y2[kEEmcNumStrips]; 
00101   float length[kEEmcNumStrips];
00102   float x0Edge[kEEmcNumEdgeStrips];
00103   float y1Edge[kEEmcNumEdgeStrips];
00104   float y2Edge[kEEmcNumEdgeStrips];
00105   float lengthEdge[kEEmcNumEdgeStrips];
00106 
00107 // fill variable arrays with data in EmcStripGeom.h 
00108   for (int i = 0; i < kEEmcNumStrips; i++) {
00109        x0[i] = EEmcStripGeomData[i].x0;
00110        y1[i] = EEmcStripGeomData[i].y1;
00111        y2[i] = EEmcStripGeomData[i].y2;
00112        length[i] = EEmcStripGeomData[i].length;
00113   }
00114   for (int i = 0; i < kEEmcNumEdgeStrips; i++) {
00115        x0Edge[i] = EEmcEdgeStripGeomData[i].x0;
00116        y1Edge[i] = EEmcEdgeStripGeomData[i].y1;
00117        y2Edge[i] = EEmcEdgeStripGeomData[i].y2;
00118        lengthEdge[i] = EEmcEdgeStripGeomData[i].length;
00119   }
00120 /*
00121 // write x, y, and length to a file for a check
00122   FILE *fp;
00123   fp = fopen("arrCheck.lis", "w");
00124   fprintf(fp, " Reqular\n");
00125   for (int i = 0; i < kEEmcNumStrips; i++) {
00126        fprintf(fp, " strip %d: x0, y1, y2, length = %f %f %f %f\n",
00127                                 i, x0[i], y1[i],  y2[i], length[i]);
00128   }
00129   fprintf(fp, " Edge\n");
00130   for (int i = 0; i < kEEmcNumEdgeStrips; i++) {
00131        fprintf(fp, " strip %d: x0, y1, y2, length = %f %f %f %f\n",
00132                   i, x0Edge[i], y1Edge[i],  y2Edge[i], lengthEdge[i]);
00133   }
00134   fclose(fp);
00135 */
00136   float delPhi = 2*pi/degree/kEEmcNumSectors;
00137   float PhiRotation[kEEmcNumSmdUVs][kEEmcNumSectors];
00138 
00139 // calculate rotation angles to prepare for local-global transformation
00140   for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
00141     for(int iSec=0; iSec<kEEmcNumSectors; iSec++){
00142        PhiRotation[iUV][iSec]=(-15.0 + iSec*delPhi)*degree; 
00143        if(iUV == 1) PhiRotation[iUV][iSec] = -1.0*PhiRotation[iUV][iSec];  
00144     }
00145   }
00146 
00147 // loop over planes
00148   for (int iPlane = 0; iPlane < kEEmcNumSmdPlanes; iPlane++) {
00149     float globalX1, globalY1, globalX2, globalY2;
00150     float x0Corr, y1Corr, y2Corr, lengthCorr; 
00151     float phi1, phi2, phiMin, phiMax;
00152     float r, rMin, rMax;
00153   
00154     mEEmcSmdParam.zPlane[iPlane] = kEEmcZSMD + 
00155                      (iPlane - kEEmcNumSmdPlanes + 2) * kEEmcSmdZPlaneShift ;
00156 // loop over UV
00157     for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) { 
00158 
00159 // loop over sectors
00160       for(int iUVSec=iPlane+1-iUV; iUVSec<kEEmcNumSectors+1-iUV; 
00161                                      iUVSec=iUVSec+kEEmcNumSmdPlanes) {
00162         int iSec;
00163         if(iUVSec == 12) iSec = 0;          
00164         else iSec = iUVSec;
00165 
00166         if(IsSectorIn(iSec)) {
00167 
00168           rMin = 1000.0;
00169           rMax = 0.0;
00170           phiMin = pi; 
00171           phiMax = -pi;
00172 
00173 // loop over strips     
00174           for (int iStrip = 0; iStrip < kEEmcNumStrips; iStrip++) {   
00175             if(kEEmcSmdMapEdge[iPlane][iSec] && iStrip > kEEmcNumEdgeStrips-1) 
00176                                                                     break;      
00177 // Id = index + 1 in all cases
00178             StructEEmcStripId  stripStructId;
00179             stripStructId.sectorId = iSec+1;
00180             stripStructId.UVId = iUV+1;
00181             stripStructId.stripId = iStrip + 1;
00182             stripStructId.planeId = iPlane + 1;
00183 
00184             StructEEmcStrip*    stripPtr = new StructEEmcStrip;
00185             stripPtr->stripStructId = stripStructId; 
00186 
00187 // correct for radius offset 
00188             x0Corr = x0[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
00189             y2Corr = y2[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
00190             if(kEEmcSmdMapEdge[iPlane][iSec]) {      
00191               y1Corr = y1Edge[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
00192               lengthCorr = lengthEdge[iStrip];
00193             }
00194             else { 
00195               y1Corr = y1[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
00196               lengthCorr = length[iStrip];
00197             }
00198 // Transform local to gloabal by rotation. Local X axis is perpendicular to
00199 // strip, Y is parallel. After rotation, x & y should be swapped for V sector. 
00200             
00201             if(iUV == 0) {
00202               globalX1 = x0Corr*cos(PhiRotation[iUV][iSec])+ 
00203                                      y1Corr*sin(PhiRotation[iUV][iSec]);
00204               globalY1 = y1Corr*cos(PhiRotation[iUV][iSec]) - 
00205                                      x0Corr*sin(PhiRotation[iUV][iSec]);
00206               globalX2 = x0Corr*cos(PhiRotation[iUV][iSec]) + 
00207                                      y2Corr*sin(PhiRotation[iUV][iSec]);
00208               globalY2 = y2Corr*cos(PhiRotation[iUV][iSec]) - 
00209                                      x0Corr*sin(PhiRotation[iUV][iSec]);
00210             }
00211             else {
00212               globalX1 = y1Corr*cos(PhiRotation[iUV][iSec]) - 
00213                                       x0Corr*sin(PhiRotation[iUV][iSec]);
00214               globalY1 = x0Corr*cos(PhiRotation[iUV][iSec]) + 
00215                                       y1Corr*sin(PhiRotation[iUV][iSec]);
00216               globalX2 = y2Corr*cos(PhiRotation[iUV][iSec]) - 
00217                                       x0Corr*sin(PhiRotation[iUV][iSec]);
00218               globalY2 = x0Corr*cos(PhiRotation[iUV][iSec]) + 
00219                                       y2Corr*sin(PhiRotation[iUV][iSec]);
00220             }
00221             r = ::sqrt(globalX1*globalX1 + globalY1*globalY1);
00222             if(r < rMin) rMin = r;
00223             r = ::sqrt(globalX2*globalX2 + globalY2*globalY2);
00224             if(r > rMax) rMax = r;
00225 
00226 //Fill StripPtrVec 
00227             stripPtr->end1.SetX(globalX1) ;
00228             stripPtr->end1.SetY(globalY1) ;
00229             stripPtr->end1.SetZ(mEEmcSmdParam.zPlane[iPlane]);
00230             stripPtr->end2.SetX(globalX2) ;
00231             stripPtr->end2.SetY(globalY2) ;
00232             stripPtr->end2.SetZ(mEEmcSmdParam.zPlane[iPlane]);
00233             stripPtr->length = lengthCorr;
00234 
00235             phi1 = stripPtr->end1.Phi();
00236             phi2 = stripPtr->end2.Phi();
00237  
00238             if(iSec !=  kEEmcSmdSectorIdPhiCrossPi - 1) {
00239               if(phi1 < phiMin) phiMin = phi1; 
00240               if(phi1 > phiMax) phiMax = phi1; 
00241               if(phi2 < phiMin) phiMin = phi2; 
00242               if(phi2 > phiMax) phiMax = phi2; 
00243             }
00244             else {
00245               if(phi1 > 0)  if(phi1 < phiMin) phiMin = phi1; 
00246               if(phi1 < 0 ) if(phi1 > phiMax) phiMax = phi1; 
00247               if(phi2 > 0)  if(phi2 < phiMin) phiMin = phi2; 
00248               if(phi2 < 0)  if(phi2 > phiMax) phiMax = phi2; 
00249             }
00250             mEEmcSector[iUV][iSec].stripPtrVec.push_back(stripPtr);
00251           } // loop over iStrip
00252 
00253 // Fill mEEmcSectors 
00254           mEEmcSector[iUV][iSec].sectorId = iSec+1;
00255           mEEmcSector[iUV][iSec].planeId = iPlane+1;
00256           mEEmcSector[iUV][iSec].phiMin = phiMin;
00257           mEEmcSector[iUV][iSec].phiMax = phiMax;
00258           mEEmcSector[iUV][iSec].rMin = rMin;
00259           mEEmcSector[iUV][iSec].rMax = rMax;
00260 
00261         } // if sector selected
00262       } // end of iUVSec loop 
00263     } // end of iUV loop
00264   } // end of iPlane loop
00265 
00266 
00267   // build revers mapping (iUV,iSec) --> iPlane
00268   memset(kEEmcSmdMap_iPlane,0,sizeof(kEEmcSmdMap_iPlane));
00269   for (int iPlane = 0; iPlane < kEEmcNumSmdPlanes; iPlane++) 
00270     for(int iSec=0; iSec<kEEmcNumSectors; iSec++) {
00271       int iuv=kEEmcSmdMapUV[iPlane][iSec];
00272       if(iuv<0 ) continue;
00273       assert(iuv<kEEmcNumSmdUVs);
00274       kEEmcSmdMap_iPlane[iuv][iSec]=iPlane;
00275     }
00276 
00277   buildStripPtrVector();
00278 
00279 } // end of buildSmdGeom 
00280 //===========================================================
00281 //===========================================================
00282 
00283 // build mStripPtrVector with getEEmcSector()
00284 void EEmcSmdGeom::buildStripPtrVector() {
00285    StructEEmcStrip *dummyStripPtr = new StructEEmcStrip;
00286    *dummyStripPtr = initStrip();
00287    EEmcStripPtrVec stripPtrVec;
00288    EEmcStripPtrVecIter p;  
00289    for(int iSec = 0; iSec < kEEmcNumSectors; iSec++) {
00290       if(mIsSectorIn[iSec]) {
00291         for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
00292           stripPtrVec = getEEmcSector(iUV,iSec).stripPtrVec;
00293           p = stripPtrVec.begin();
00294           int PlaneId = getEEmcSector(iUV,iSec).planeId;
00295           while(p !=stripPtrVec.end()) {
00296              mStripPtrVector.push_back(*p);
00297              p++;
00298 
00299           } 
00300           if(kEEmcSmdMapEdge[PlaneId-1][iSec]) {
00301                for(int i=0; i < kEEmcNumStrips - kEEmcNumEdgeStrips; i++)   
00302                   mStripPtrVector.push_back(dummyStripPtr);
00303           } 
00304         }
00305       }
00306       else {
00307         for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
00308           for(int iStrip = 0; iStrip < kEEmcNumStrips; iStrip++) 
00309              mStripPtrVector.push_back(dummyStripPtr);
00310         }
00311       }
00312    }
00313 }
00314 
00315 // set status of selected sectors
00316 void EEmcSmdGeom::setSectors(const intVec sectorIdVec) {
00317        for(int iSec = 0; iSec< kEEmcNumSectors; iSec++) 
00318                                    mIsSectorIn[iSec] = false;
00319        for(unsigned int i = 0; i < sectorIdVec.size(); i++) { 
00320           for(int iSec = 0; iSec< kEEmcNumSectors; iSec++) {
00321             if (sectorIdVec[i] == iSec+1) mIsSectorIn[iSec] = true;
00322           }
00323        }                                       
00324 }
00325 
00326 // instance and initialize a strip
00327 StructEEmcStrip EEmcSmdGeom::initStrip() const {
00328     TVector3  zero = 0.0;
00329     StructEEmcStrip strip; 
00330     strip.stripStructId.stripId = 0;
00331     strip.stripStructId.UVId = 0;
00332     strip.stripStructId.sectorId = 0;
00333     strip.stripStructId.planeId = 0;
00334     strip.end1 = zero;
00335     strip.end2 = zero;
00336     strip.length = 0.0;
00337     return strip;
00338 }
00339 
00340 // return index of a sector from a global point in a plane 
00341 Int_t EEmcSmdGeom::getEEmcISec(const Int_t iPlane, 
00342                            const TVector3& point) const {
00343      Int_t indexSec = -1;
00344      float phiMin, phiMax, rMin, rMax;
00345      float phi = point.Phi();
00346      float r = ::sqrt(point.x()*point.x() + point.y()*point.y());
00347 
00348      for (int iSec = 0; iSec < kEEmcNumSectors; iSec++) {
00349        int iUV = kEEmcSmdMapUV[iPlane][iSec];
00350        if(iUV >= 0 && IsSectorIn(iSec)) {       
00351          phiMin = mEEmcSector[iUV][iSec].phiMin;  
00352          phiMax = mEEmcSector[iUV][iSec].phiMax;  
00353          rMin = mEEmcSector[iUV][iSec].rMin;  
00354          rMax = mEEmcSector[iUV][iSec].rMax;
00355          if(iSec !=  kEEmcSmdSectorIdPhiCrossPi - 1) {
00356            if (phi >= phiMin && phi < phiMax && r > rMin && r < rMax) {
00357               indexSec = iSec;
00358               break;
00359            }
00360          }
00361 // sector9 between 165 deg (Min) and -165 deg (Max)
00362          else { 
00363            if(((phi > 0.0 && phi >= phiMin) || (phi < 0.0 && phi < phiMax))
00364                                                  && r > rMin && r < rMax){
00365                indexSec = iSec;
00366                break;
00367            }
00368          }
00369        }
00370      }
00371      return indexSec; 
00372 }
00373 
00374 // return a strip pointer from indices   
00375 StructEEmcStrip* EEmcSmdGeom::getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec) {
00376     int i = iStrip + iUV*kEEmcNumStrips + iSec*kEEmcNumStrips*kEEmcNumSmdUVs;
00377     return mStripPtrVector[i];
00378 }
00379 
00380 // return a strip pointer from indices   
00381 const StructEEmcStrip* EEmcSmdGeom::getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec) const {
00382     int i = iStrip + iUV*kEEmcNumStrips + iSec*kEEmcNumStrips*kEEmcNumSmdUVs;
00383     return mStripPtrVector[i];
00384 }
00385 
00386 
00387 //==================================================================
00388 // get DCA strip pointer from a point  
00389 // iPlane=[0,1,2] - experts only, changes meaning form sector to sector
00390 const StructEEmcStrip* EEmcSmdGeom::getDcaStripPtr(const Int_t iPlane, 
00391                const Int_t iSec, const TVector3& point, Float_t* dca) const
00392 {
00393   //    StructEEmcStrip* stripPtr;
00394   //    stripPtr = new StructEEmcStrip;
00395     *dca = 1000.0;
00396     int iStrip = -1;
00397     //$$$    int iUV;
00398     float x1,y1,x2,y2,mu,d;
00399     EEmcStripPtrVec stripPtrVec;  
00400     EEmcStripPtrVecIter p;  
00401 
00402     
00403     Int_t iUV = kEEmcSmdMapUV[iPlane][iSec]; // jcw 2/6/04 moved here from $$$
00404                                              // iUV may be used uninitialized
00405                                              // below, otherwise...
00406 //  int iSec = getEEmcISec(iPlane, point);
00407     if(iSec >= 0 && IsSectorIn(iSec)) {
00408       //$$$iUV = kEEmcSmdMapUV[iPlane][iSec];
00409       stripPtrVec =  getEEmcSector(iUV,iSec).stripPtrVec;
00410       p =  stripPtrVec.begin();
00411       while(p != stripPtrVec.end()) {
00412         x1 = (*p)->end1.x();
00413         y1 = (*p)->end1.y();
00414         x2 = (*p)->end2.x();
00415         y2 = (*p)->end2.y();
00416         mu = -1.0/::sqrt((y2-y1)*(y2-y1) + (x2-x1)*(x2-x1)) *
00417           ((x2*y1-x1*y2)/fabs(x2*y1-x1*y2));
00418 // distance d carries a sign
00419         d = ((y2-y1)*point.x() + (x1-x2)*point.y() + (x2*y1-x1*y2))*mu;
00420         
00421         if(fabs(d) < fabs(*dca)) {
00422           *dca = d;
00423           iStrip = (*p)->stripStructId.stripId - 1;
00424         }
00425         if(d < 0) break;
00426         p++;
00427       }
00428     }
00429     if(iStrip >=0) {
00430       //stripPtr = getStripPtr(iStrip,iUV,iSec);
00431       //return stripPtr;
00432       return getStripPtr(iStrip,iUV,iSec);
00433     }
00434     else {
00435       StructEEmcStrip *stripPtr = new StructEEmcStrip;
00436       (*stripPtr) = initStrip();
00437       //    std::cout << "NO dca strip found in plane (sector empty or not in)" << std::endl;  //silentium
00438 
00439       return stripPtr;
00440     }
00441 }
00442 
00443 //==================================================================
00444 // get DCA strip pointer from a point  
00445 // iPlane=[0,1,2] - experts only, changes meaning form sector to sector
00446 const StructEEmcStrip* EEmcSmdGeom::getDcaStripPtr(const Int_t iPlane, 
00447                          const TVector3& point, Float_t* dca) const {
00448   int iSec = getEEmcISec(iPlane, point);
00449   return  getDcaStripPtr(iPlane, iSec, point, dca); 
00450 }
00451 
00452 
00453 //==================================================================
00454 // get DCA strip pointer from a point  
00455 // iUV=[0,1] , maps [U,V]
00456 // Warn, this code does not handle well sector boundaries, it ignores possible tripple overlaps and uses only nominal U or V planes in a given sector
00457 const StructEEmcStrip* EEmcSmdGeom::getDca2Strip(const Int_t iUV, 
00458                                            const TVector3& point, Float_t* dca) const {
00459   assert(iUV>=0 || iUV<kEEmcNumSmdUVs); 
00460   float phiDeg=atan2(point.y(),point.x())/3.1316*180.;
00461   //printf("phiDeg=%.1f  \n",phiDeg);
00462   int iSec= ((int) ( 12.-(phiDeg-75.)/30.) )%12;
00463   assert(iSec>=0);
00464   assert( iSec<kEEmcNumSectors);
00465   int iPlane= kEEmcSmdMap_iPlane[iUV][iSec];// now find  mapping iUV --> iPlane
00466   assert(iPlane>=0 && iPlane<kEEmcNumSmdPlanes);
00467   return  getDcaStripPtr(iPlane, iSec, point, dca); 
00468 }
00469 
00470 
00471 //==================================================================
00472 // match two strips  
00473   bool EEmcSmdGeom::matchStrips(const StructEEmcStripId &stripStructId1, 
00474                                       const StructEEmcStripId &stripStructId2,
00475                                       Int_t nTolerance) const {
00476     bool match = false;
00477     if(stripStructId1.UVId == stripStructId2.UVId &&
00478        stripStructId1.sectorId == stripStructId2.sectorId) {
00479          if((TMath::Abs(stripStructId1.stripId - stripStructId2.stripId) <= nTolerance))
00480            match = true;
00481     }
00482     return match;
00483 }
00484 
00485 
00486 
00487 TVector3  EEmcSmdGeom::getstripEnd(const StructEEmcStrip &strip, 
00488                                                     const Int_t endId) const {
00489       TVector3 end;
00490       if(endId == 1) end = strip.end1;
00491       else end = strip.end2;
00492 
00493       return end;
00494 }
00495 
00496 // methods of printout
00498 void EEmcSmdGeom::printGeom(ostream& os) const {
00499   os << "------EEmcSmdGeom::printGeom()------" << endl;
00500   os << " " << "z[3]          = " 
00501      << " " << getEEmcSmdParam().zPlane[0] 
00502      << " " << getEEmcSmdParam().zPlane[1] 
00503      << " " << getEEmcSmdParam().zPlane[2] << endl;
00504   os << " " << "rOffset[3]    = "
00505      << " " << getEEmcSmdParam().rOffset[0]
00506      << " " << getEEmcSmdParam().rOffset[1]
00507      << " " << getEEmcSmdParam().rOffset[2] << endl;
00508   os << " " << "stripWidth    = "
00509      << " " << getEEmcSmdParam().stripWidth << endl;
00510   os << "---------------------------------------" << endl;
00511 }
00512 
00514 void EEmcSmdGeom::printSector(const StructEEmcSmdSector sector, ostream& os) const {
00515   float delPhi;
00516   int iUV = kEEmcSmdMapUV[sector.planeId-1][sector.sectorId-1];
00517   delPhi = (sector.phiMax - sector.phiMin)/degree;
00518   if(sector.sectorId == kEEmcSmdSectorIdPhiCrossPi) 
00519                              delPhi = 2*pi/degree + delPhi;
00520   
00521   os << "------EEmcSmdGeom::printSector()------" << endl;
00522   os << kEEmcSmdUVChar[iUV] << " Sector:  sectorId, planeId, nStrips      = " 
00523      << " " << sector.sectorId 
00524      << " " << sector.planeId 
00525      << " " << sector.stripPtrVec.size() << endl;
00526   os << "           phiMin, phiMax, delPhi  = " 
00527      << " " << sector.phiMin/degree 
00528      << " " << sector.phiMax/degree
00529      << " " << delPhi << endl;
00530   os << "           rMin, rMax delR         = "  
00531      << " " << sector.rMin 
00532      << " " << sector.rMax 
00533      << " " << sector.rMax - sector.rMin << endl;
00534   os << "------------------------------------" << endl;
00535 }
00536 
00538 void EEmcSmdGeom::printStrip(const StructEEmcStrip strip, ostream& os) const {
00539   char UVChar;  
00540   if(strip.stripStructId.sectorId == 0) UVChar = 'X';
00541   else
00542     UVChar = kEEmcSmdUVChar[strip.stripStructId.UVId - 1];
00543 
00544   os << "------EEmcSmdGeom::printStrip()------" << endl;
00545 
00546     os << "Strip:  sectorId, planeUV, stripId, planeId    = "
00547        << " " << strip.stripStructId.sectorId
00548        << " " << UVChar 
00549        << " " << strip.stripStructId.stripId
00550        << " " << strip.stripStructId.planeId << endl;
00551     os << "        x1, y1, x2, y2, z     = "
00552        << " " << strip.end1.x() 
00553        << " " << strip.end1.y() 
00554        << " " << strip.end2.x() 
00555        << " " << strip.end2.y() 
00556        << " " << strip.end2.z() << endl; 
00557     os << "        phi1, phi2, length    = "
00558        << " " << strip.end1.Phi()/degree
00559        << " " << strip.end2.Phi()/degree
00560        << " " << strip.length << endl;
00561     os << "------------------------------------" << endl;
00562 }
00563 
00565 void EEmcSmdGeom::printStripId(const StructEEmcStripId stripStructId, ostream& os) const {
00566   char UVChar;  
00567   if(stripStructId.sectorId == 0) UVChar = 'X';
00568   else
00569     UVChar = kEEmcSmdUVChar[stripStructId.UVId - 1];
00570 
00571   os << "------EEmcSmdGeom::printStripId()------" << endl;
00572     os << "Strip:  sectorId, stripId, planeId    = "
00573        << " " << stripStructId.sectorId
00574               << UVChar 
00575        << " " << stripStructId.stripId
00576        << " " << stripStructId.planeId << endl;
00577     os << "------------------------------------" << endl;
00578 }
00579 
00580 
00581 
00582 // Move to StEEmcSmdGeom
00583 #if 0
00584 // printout delPhi and centerPhi used in ITTF
00585 void EEmcSmdGeom::printSectorPhis(const Int_t iPlane, const Int_t iSec,
00586                                                               ostream& os ) {
00587   int iUV;
00588   iUV = kEEmcSmdMapUV[iPlane][iSec];
00589 
00590   os << "------EEmcSmdGeom::printPhis()------" << endl;
00591   os << " planeId = " << iPlane + 1 << " sectorId = " << iSec + 1 << endl;
00592   if(iUV >= 0) 
00593     os << " " <<  kEEmcSmdUVChar[iUV] << " Sector" << endl; 
00594   else  
00595     os << " Empty" << endl; 
00596   os << " delPhi = " << getEEmcSmdDelPhi(iPlane, iSec)/degree <<
00597     " " << "centerPhi = " << getEEmcSmdCenterPhi(iPlane, iSec)/degree 
00598      << endl;
00599   
00600 }
00601 #endif
00602 
00603 
00605 
00606 //
00607 // Function(s) to return the position of the crossing of two strips, given
00608 //   either the sectorID and strip ID's, or the StructEEmcStrips.
00609 //
00610 
00611 // Wrapper function when we don't have the actual strips, just sector and
00612 //   strip ID's.
00613 TVector3 EEmcSmdGeom::getIntersection ( Int_t sector, Int_t uId, Int_t vId ) const {
00614 
00615   Int_t nU = getNStrips( sector, 0 );
00616   Int_t nV = getNStrips( sector, 1 );
00617   if ( uId >= nU || vId >= nV ) {
00618     LOG_DEBUG<<"::getIntersection(...) passed invalid strip ID sector="<<sector<<" uId="<<uId<<" vId="<<vId<<std::endl;
00619     return TVector3(1.,1.,-999.0);
00620   }
00621 
00622   return getIntersection ( getStripPtr( uId, 0, sector ),
00623                            getStripPtr( vId, 1, sector ) );
00624 
00625 }
00626 
00627 
00628 TVector3 EEmcSmdGeom::getIntersection ( Int_t sector, Int_t uId, Int_t vId, const TVector3 &vert ) const {
00629 
00630   Int_t nU = getNStrips( sector, 0 );
00631   Int_t nV = getNStrips( sector, 1 );
00632   if ( uId >= nU || vId >= nV ) {
00633     LOG_DEBUG<<"::getIntersection(...) passed invalid strip ID sector="<<sector<<" uId="<<uId<<" vId="<<vId<<std::endl;
00634     return TVector3(1.,1.,-999.0);
00635   }
00636 
00637   return getIntersection ( getStripPtr( uId, 0, sector ),
00638                            getStripPtr( vId, 1, sector ), vert );
00639 
00640 }
00641 
00642 
00643 /*********************************************************/ 
00644 TVector3 EEmcSmdGeom::getIntersection ( const StructEEmcStrip *u,
00645                                         const StructEEmcStrip *v,
00646                                         const TVector3 &vertex ) const
00647 {
00648   // The strips are arranged sensibly, so that the ID's and the
00649   // widths of the strips basically tell us the position of the
00650   // crossing point.  However, we need to know a few things:
00651                                                                              
00652   Int_t uSectorId = u -> stripStructId.sectorId;   // This would be easier if
00653   Int_t vSectorId = v -> stripStructId.sectorId;   // we were dealing with a
00654   //Int_t uPlaneId  = u -> stripStructId.planeId;  // class instead of a struct
00655   //Int_t vPlaneId  = v -> stripStructId.planeId;  //
00656 
00657   Int_t uId = u -> stripStructId.stripId;
00658   Int_t vId = v -> stripStructId.stripId;
00659 
00660   // Get vectors pointing to the start of the u and v strips in question,
00661   //   as well as the ends.  NOTE: We pass uId - 1 to getStripPtr, because
00662   //   that routine expects the c++ _index_... (the convention in this
00663   //   class is that quantities beginning with an "i" correspond to a
00664   //   c++ index numbered from 0, rather than a fortran index numbered
00665   //   from 1...)
00666   TVector3 u0 = getStripPtr ( uId-1, 0, uSectorId-1 ) -> end1;
00667   TVector3 v0 = getStripPtr ( vId-1, 1, vSectorId-1 ) -> end1;
00668   
00669 
00670   TVector3 uF = getStripPtr ( uId-1, 0, uSectorId-1 ) -> end2;
00671   TVector3 vF = getStripPtr ( vId-1, 1, vSectorId-1 ) -> end2;
00672   
00673 
00674   TVector3 u0p=(u0-vertex);  // uOp is vector from vertex to beginning of first strip from vertex 
00675   TVector3 uFp=(uF-vertex);  // uFp is vector from vertex to end of first strip from vertex
00676   TVector3 v0p=(v0-vertex);  // vOp is vector from vertex to beginning of first strip from vertex
00677   TVector3 vFp=(vF-vertex);  // vFp is vector from vertex to end of first strip from vertex
00678 
00680   //          U then V!  return point on V strip
00682 
00683   // N is normal vector to plane
00684   TVector3 Nv = u0p.Cross(uFp);
00685 
00686   // use components of N to establish D in plane of form A(x-x0)+B(y-y0)+C(z-z0)+D=0, use vertex for point
00687   Float_t Dv = (Nv.X()*vertex.X() + Nv.Y()*vertex.Y() + Nv.Z()*vertex.Z());
00688 
00689   // scale determines how far along vector from v0 to vF to go to get point 
00690   // of intersection with plane.  If it's <0 or >1, then the point of 
00691   // intersection is off the physical size of the strip.  Return an error
00692   // value
00693   Float_t scale_numerv = (Nv.X()*v0.X() + Nv.Y()*v0.Y() + Nv.Z()*v0.Z() - Dv);
00694   Float_t scale_denomv = (Nv.X()*(v0.X()-vF.X()) + Nv.Y()*(v0.Y()-vF.Y()) + Nv.Z()*(v0.Z()-vF.Z()));
00695   Float_t scalev=scale_numerv/scale_denomv;
00696   if (scalev < 0 || scalev > 1) {
00697     TVector3 ErrorVector(1.,1.,-999.);
00698     LOG_DEBUG<<GetName()<<"::getIntersection( ) passed non-intersecting SMD strips " << *u << " " << *v << endm;
00699     return ErrorVector;
00700   }  
00701   
00702   
00703   // Rv is the vector to the final point from the origin ON V STRIP
00704   TVector3 Rv = (v0 + scalev*(vF-v0));
00705 
00706 
00707     
00708   // S is the vector that would go from vertex to end point
00709   // TVector3 S = (R-z0);
00710   
00712   //       V THEN U! returns point on u strip
00714   // N is normal vector to plane
00715   TVector3 Nu = v0p.Cross(vFp);
00716   
00717   // use components of N to establish D in plane of form A(x-x0)+B(y-y0)+C(z-z0)+D=0, use vertex for point
00718   Float_t Du = (Nu.X()*vertex.X() + Nu.Y()*vertex.Y() + Nu.Z()*vertex.Z());
00719   
00720   // scale determines how far along vector from v0 to vF to go to get point of intersection with plane
00721   Float_t scale_numeru = (Nu.X()*u0.X() + Nu.Y()*u0.Y() + Nu.Z()*u0.Z() - Du);
00722   Float_t scale_denomu = (Nu.X()*(u0.X()-uF.X()) + Nu.Y()*(u0.Y()-uF.Y()) + Nu.Z()*(u0.Z()-uF.Z()));
00723   Float_t scaleu=scale_numeru/scale_denomu;
00724   if (scaleu < 0 || scaleu > 1) {
00725     TVector3 ErrorVector(1.,1.,-999.);
00726     return ErrorVector;
00727   }
00728   
00729   
00730   // Ru is the vector to the final point from the origin ON U STRIP
00731   TVector3 Ru = (u0 + scaleu*(uF-u0));
00732   
00733   
00734   // gives us vector direction FROM THE VERTEX
00735   TVector3 R = ((Rv+Ru)*0.5);
00736   
00737   return R;
00738                                                                                                                                                   
00739 
00740 
00741 }
00742 
00743 /*********************************************************/ 
00744 
00745 
00746 TVector3 EEmcSmdGeom::getIntersection ( const StructEEmcStrip *u,
00747                                         const StructEEmcStrip *v
00748                                         ) const {
00749 
00750   return getIntersection( u, v, TVector3(0,0,0) );
00751 
00752 }
00753 
00754 
00755 /****************************************************************************/
00756 
00757 // output strip for smd strips
00758 ostream& operator<<(ostream &os, const StructEEmcStrip &strip)
00759 {
00760   const Char_t *nameUV[]={"U","V"};
00761   TString stripname="";
00762   if ( strip.stripStructId.sectorId < 10 ) stripname+="0"; stripname+=strip.stripStructId.sectorId;
00763   stripname += nameUV[ strip.stripStructId.UVId-1 ];
00764   if ( strip.stripStructId.stripId < 10 ) stripname+="0";
00765   if ( strip.stripStructId.stripId < 100 ) stripname+="0";
00766   stripname+=strip.stripStructId.stripId;
00767   os << stripname << " depth=" << strip.stripStructId.planeId;
00768   return os;
00769 }
00770 
00771 
00773 /*
00774  * $Log: EEmcSmdGeom.cxx,v $
00775  * Revision 1.16  2010/08/26 22:48:55  ogrebeny
00776  * Improved constness
00777  *
00778  * Revision 1.15  2009/08/25 18:33:11  fine
00779  * fix the compilation issues under SL5_64_bits  gcc 4.3.2
00780  *
00781  * Revision 1.14  2008/02/17 17:37:32  balewski
00782  * demote warning about strips intersecting beyond modle to DEBUG
00783  *
00784  * Revision 1.13  2007/07/12 19:30:15  fisyak
00785  * Add includes for ROOT 5.16
00786  *
00787  * Revision 1.12  2007/05/01 20:22:25  jwebb
00788  * Added error handling to EEmcSmdGeom::getIntersection().
00789  *
00790  * Revision 1.11  2007/02/17 01:29:05  balewski
00791  * less printout
00792  *
00793  * Revision 1.10  2007/02/02 02:11:11  balewski
00794  * simplification of  EEmcSmdGeom::getDca2Strip(..)
00795  *
00796  * Revision 1.9  2007/02/01 20:33:46  wzhang
00797  * dca initialized in getDcaStripPtr()
00798  *
00799  * Revision 1.8  2007/02/01 13:47:39  balewski
00800  * bug fix in getDca2Strip(), more methodhs are public
00801  *
00802  * Revision 1.7  2007/01/26 00:51:08  balewski
00803  * too strong protection
00804  *
00805  * Revision 1.6  2007/01/25 22:33:21  balewski
00806  * add:
00807  * - better writeup
00808  * - new simpler to use method calculating dca fo track to strip, it is just a wrapper, some approximations were used, may fail at the sector boundary
00809  *
00810  * Revision 1.5  2007/01/12 23:53:14  jwebb
00811  * Fix applied to eliminate parralax error in the EEmcSmdGeom::getIntersection()
00812  * method.
00813  *
00814  * Revision 1.4  2004/06/03 23:01:06  jwebb
00815  * Fixed memory leak reported by Bob.  Fixed another memory leak.
00816  *
00817  * Revision 1.3  2004/02/06 22:33:08  jwebb
00818  * Moved statement to fix warning.
00819  *
00820  * Revision 1.2  2004/01/29 16:37:25  jwebb
00821  * Removed dependence on StMaker.h and PhysicalConstants.h.  Should be fully
00822  * decoupled from Star environment now.
00823  *
00824  * Revision 1.1  2004/01/29 15:26:10  jwebb
00825  * The StEEmcSmdGeom class was split into two classes.  All StRoot-independent
00826  * code has been moved to EEmcSmdGeom.  TVector3 replaces StThreeVectorD in
00827  * all function calls in EEmcSmdGeom.  StThreeVectorD wrappers are provided
00828  * in StEEmcSmdGeom, for integration into Star framework.
00829  *
00830  *
00831  *****************************************************************************
00832  * Log: StEEmcSmdGeom.cxx,v 
00833  * Revision 1.8  2004/01/12 14:34:09  wzhang
00834  * Corrected the usage of EEmcStripPtrVecIter
00835  *
00836  * Revision 1.7  2003/12/05 00:06:10  jwebb
00837  * Member function added to return a vector pointing to the intersection of
00838  * two strips.
00839  *
00840  * Revision 1.6  2003/10/15 15:26:08  wzhang
00841  * improved and reorganized
00842  *
00843  * Revision 1.5  2003/09/02 17:57:56  perev
00844  * gcc 3.2 updates + WarnOff
00845  *
00846  * Revision 1.4  2003/08/22 15:14:26  wzhang
00847  * Added ClassImp and method stripEnd
00848  *
00849  * Revision 1.3  2003/06/11 18:58:19  wzhang
00850  * added geometry methods for StiEEmc
00851  *
00852  * Revision 1.2  2003/04/04 15:33:32  wzhang
00853  * included EEmcGeomDefs.h & improved codes
00854  *
00855  * Revision 1.1  2003/03/28 15:49:59  balewski
00856  * first
00857  * 
00858  *
00859  *******************************************************************/