00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <iostream>
00013 #include <vector>
00014 #include "TGeoBBox.h"
00015
00016 #include "MucGeomSvc/MucGeoGap.h"
00017 #include "MucGeomSvc/MucGeoStrip.h"
00018 #include "MucGeomSvc/MucGeometron.h"
00019 #include "MucGeomSvc/MucConstant.h"
00020 #include "Identifier/MucID.h"
00021
00022 using namespace std;
00023
00024 MucGeoGap::MucGeoGap()
00025 : m_Part(-1),
00026 m_Seg(-1),
00027 m_Gap(-1),
00028 m_StripNum(-1),
00029 m_XSize(-1.0),
00030 m_YSize(-1.0),
00031 m_ZSize(-1.0)
00032 {
00033
00034 }
00035
00036 MucGeoGap::MucGeoGap(const int part,
00037 const int seg,
00038 const int gap,
00039 const int orient,
00040 const int stripNum,
00041 const float xSize,
00042 const float ySize,
00043 const float zSize,
00044 const float xTarget1Global,
00045 const float yTarget1Global,
00046 const float zTarget1Global,
00047 const float xTarget2Global,
00048 const float yTarget2Global,
00049 const float zTarget2Global,
00050 const float xTarget3Global,
00051 const float yTarget3Global,
00052 const float zTarget3Global,
00053 const float dzHighEdge,
00054 const float dzFarFrontGas,
00055 const float dzNearFrontGas,
00056 const float dzNearBackGas,
00057 const float dzFarBackGas,
00058 const float dxTarget1ToFiducial,
00059 const float dyTarget1ToFiducial,
00060 const float dxFiducialToCenter,
00061 const float dyFiducialToCenter)
00062 : m_Part(part),
00063 m_Seg(seg),
00064 m_Gap(gap),
00065 m_StripNum(stripNum),
00066 m_Orient(orient),
00067 m_XSize(xSize),
00068 m_YSize(ySize),
00069 m_ZSize(zSize),
00070 m_dzHighEdge(dzHighEdge),
00071 m_dzFarFrontGas(dzFarFrontGas),
00072 m_dzNearFrontGas(dzNearFrontGas),
00073 m_dzNearBackGas(dzNearBackGas),
00074 m_dzFarBackGas(dzFarBackGas),
00075 m_pMucGeoStrip(MucID::getStripNum(part, seg, gap))
00076 {
00077
00078
00079
00080
00081
00082 Hep3Vector v1Global(xTarget1Global, yTarget1Global, zTarget1Global);
00083 Hep3Vector v2Global(xTarget2Global, yTarget2Global, zTarget2Global);
00084 Hep3Vector v3Global(xTarget3Global, yTarget3Global, zTarget3Global);
00085
00086 Hep3Vector v2To1Global = v1Global - v2Global;
00087 Hep3Vector v3To1Global = v1Global - v3Global;
00088
00089
00090 Hep3Vector v1Gap(dxTarget1ToFiducial + dxFiducialToCenter,
00091 dyTarget1ToFiducial + dyFiducialToCenter,
00092 0.0);
00093
00094
00095 Hep3Vector newX(1,0,0), newY(0,1,0), newZ(0,0,1);
00096 HepRotation rotateGap(newX, newY, newZ);
00097
00098 if(part == 1) {
00099 newX = v2To1Global;
00100 newX.setMag(1.0);
00101 newY = v3To1Global;
00102 newY.setMag(1.0);
00103 newZ = newX.cross(newY);
00104 }
00105 if(part != 1) {
00106 newY = -1.0*v2To1Global;
00107 newY.setMag(1.0);
00108 newZ = v3To1Global.cross(newY);
00109 newZ.setMag(1.0);
00110 newX = newY.cross(newZ);
00111 }
00112
00113 HepRotation rotateGlobal(newX, newY, newZ);
00114 HepRotation rotateGapToGlobal(rotateGlobal * rotateGap);
00115
00116 Hep3Vector translateGapToGlobal(v1Global - rotateGapToGlobal * v1Gap);
00117
00118 m_Rotation = rotateGapToGlobal;
00119 m_Translation = translateGapToGlobal;
00120 m_Center = m_Translation;
00121
00122
00123
00124
00125 HepMatrix rM(3,3), rMT(3,3);
00126 for(int i=0; i<3; i++){
00127 for(int j=0; j<3; j++){
00128 rM[i][j] = m_Rotation(i,j);
00129 }
00130 }
00131 rMT = rM.T();
00132
00133 Hep3Vector newXT(rMT(1,1),rMT(2,1),rMT(3,1));
00134 Hep3Vector newYT(rMT(1,2),rMT(2,2),rMT(3,2));
00135 Hep3Vector newZT(rMT(1,3),rMT(2,3),rMT(3,3));
00136 HepRotation rotateGlobalToGap(newXT,newYT,newZT);
00137
00138 m_RotationT = rotateGlobalToGap;
00139 }
00140
00142 MucGeoGap::MucGeoGap(const int part,
00143 const int seg,
00144 const int gap,
00145 const int orient,
00146 const int stripNum,
00147 const TGeoPhysicalNode *gapPhysicalNode,
00148 const float ironThickness)
00149 : m_Part(part),
00150 m_Seg(seg),
00151 m_Gap(gap),
00152 m_StripNum(stripNum),
00153 m_Orient(orient),
00154 m_pMucGeoStrip(MucID::getStripNum(part, seg, gap)),
00155 m_IronThickness(ironThickness)
00156
00157 {
00158 TGeoBBox *gapBox = (TGeoBBox*)gapPhysicalNode->GetShape();
00159
00160 m_XSize = gapBox->GetDX() * 2.0;
00161 m_YSize = gapBox->GetDY() * 2.0;
00162 m_ZSize = gapBox->GetDZ() * 2.0;
00163
00164 m_dzHighEdge = 20.0;
00165 m_dzFarFrontGas = -4.5;
00166 m_dzNearFrontGas = -2.5;
00167 m_dzNearBackGas = 2.5;
00168 m_dzFarBackGas = 4.5;
00169
00170 double eRot[9], *pRot;
00171 pRot = eRot;
00172 pRot = gapPhysicalNode->GetMatrix()->GetRotationMatrix();
00173
00174
00175
00176
00177
00178 Hep3Vector rotX(pRot[0], pRot[3], pRot[6]);
00179 Hep3Vector rotY(pRot[1], pRot[4], pRot[7]);
00180 Hep3Vector rotZ(pRot[2], pRot[5], pRot[8]);
00181 HepRotation rotateGapToGlobal(rotX, rotY, rotZ);
00182 m_Rotation = rotateGapToGlobal;
00183
00184 Hep3Vector rotTX(pRot[0], pRot[1], pRot[2]);
00185 Hep3Vector rotTY(pRot[3], pRot[4], pRot[5]);
00186 Hep3Vector rotTZ(pRot[6], pRot[7], pRot[8]);
00187 HepRotation rotateGlobalToGap(rotTX, rotTY, rotTZ);
00188 m_RotationT = rotateGlobalToGap;
00189
00190 double eTrans[3], *pTrans;
00191 pTrans = eTrans;
00192 pTrans = gapPhysicalNode->GetMatrix()->GetTranslation();
00193
00194 Hep3Vector translateGapToGlobal(pTrans[0], pTrans[1], pTrans[2]);
00195 m_Translation = translateGapToGlobal;
00196 m_Center = m_Translation;
00197
00198
00199 double eTrans_gap[3], *pTrans_gap;
00200 pTrans_gap = eTrans_gap;
00201 pTrans_gap = gapPhysicalNode->GetMatrix(2)->GetTranslation();
00202
00203
00204 Hep3Vector GapSurface1, GapSurface2;
00205 Hep3Vector GapCenter(pTrans_gap[0], pTrans_gap[1], pTrans_gap[2]);
00206 GapSurface1 = GapCenter - rotZ * 20;
00207 GapSurface2 = GapCenter + rotZ * 20;
00208
00209 m_SurfaceInner = GapSurface1;
00210 m_SurfaceOuter = GapSurface2;
00211 if(GapSurface1.mag()>GapSurface2.mag())
00212 {
00213 m_SurfaceInner = GapSurface2;
00214 m_SurfaceOuter = GapSurface1;
00215 }
00216
00217
00218
00219
00220 }
00221
00222 MucGeoGap&
00223 MucGeoGap::operator=(const MucGeoGap &orig)
00224 {
00225
00226 if (this != &orig) {
00227 m_Part = orig.m_Part;
00228 m_Seg = orig.m_Seg;
00229 m_Gap = orig.m_Gap;
00230 m_StripNum = orig.m_StripNum;
00231 m_Orient = orig.m_Orient;
00232 m_HitStatus = orig.m_HitStatus;
00233
00234 m_XSize = orig.m_XSize;
00235 m_YSize = orig.m_YSize;
00236 m_ZSize = orig.m_ZSize;
00237
00238 m_dzHighEdge = orig.m_dzHighEdge;
00239 m_dzFarFrontGas = orig.m_dzFarFrontGas;
00240 m_dzNearFrontGas = orig.m_dzNearFrontGas;
00241 m_dzNearBackGas = orig.m_dzNearBackGas;
00242 m_dzFarBackGas = orig.m_dzFarBackGas;
00243
00244
00245 m_Center = orig.m_Center;
00246 m_Rotation = orig.m_Rotation;
00247 m_RotationT = orig.m_RotationT;
00248 m_Translation = orig.m_Translation;
00249
00250 m_IronThickness = orig.m_IronThickness;
00251
00252 m_pMucGeoStrip = orig.m_pMucGeoStrip;
00253 }
00254 return *this;
00255 }
00256
00257 MucGeoGap::MucGeoGap(const MucGeoGap &orig)
00258 : m_Part(orig.m_Part),
00259 m_Seg(m_Seg),
00260 m_Gap(m_Gap),
00261 m_StripNum(orig.m_StripNum),
00262 m_Orient(orig.m_Orient),
00263 m_HitStatus(orig.m_HitStatus),
00264 m_XSize(orig.m_XSize),
00265 m_YSize(orig.m_YSize),
00266 m_ZSize(orig.m_ZSize),
00267 m_dzHighEdge(orig.m_dzHighEdge),
00268 m_dzFarFrontGas(orig.m_dzFarFrontGas),
00269 m_dzNearFrontGas(orig.m_dzNearFrontGas),
00270 m_dzNearBackGas(orig.m_dzNearBackGas),
00271 m_dzFarBackGas(orig.m_dzFarBackGas),
00272 m_Center(orig.m_Center),
00273 m_Rotation(orig.m_Rotation),
00274 m_RotationT(orig.m_RotationT),
00275 m_Translation(orig.m_Translation),
00276 m_IronThickness(orig.m_IronThickness)
00277 {
00278
00279 m_pMucGeoStrip = orig.m_pMucGeoStrip;
00280 }
00281
00282 MucGeoGap::~MucGeoGap()
00283 {
00284
00285 MucGeoStrip* pStrip = 0;
00286 while(m_pMucGeoStrip.size() > 0){
00287 pStrip = m_pMucGeoStrip.back();
00288 delete pStrip;
00289 m_pMucGeoStrip.pop_back();
00290 }
00291 }
00292
00293 MucGeoStrip*
00294 MucGeoGap::GetStrip(const int strip) const
00295 {
00296
00297 if ( strip < 0 || strip >= int(m_pMucGeoStrip.size()) ) {
00298 std::cout << "Error: MucGeoGap::GetStrip() strip " << strip << " exceed strip range" << endl;
00299 return 0;
00300 }
00301 if ( !m_pMucGeoStrip[strip] ) {
00302 std::cout << "Error: MucGeoGap::GetStrip() strip " << strip << " not found" << endl;
00303 return 0;
00304 }
00305
00306 return m_pMucGeoStrip[strip];
00307 }
00308
00309 void
00310 MucGeoGap::GetSize(float &xSize, float &ySize, float &zSize) const
00311 {
00312
00313 xSize = m_XSize;
00314 ySize = m_YSize;
00315 zSize = m_ZSize;
00316
00317 return;
00318 }
00319
00320 void
00321 MucGeoGap::GetRotationMatrix(float &thetaX, float &phiX,
00322 float &thetaY, float &phiY,
00323 float &thetaZ, float &phiZ) const
00324 {
00325
00326 Hep3Vector x(m_Rotation.colX());
00327 Hep3Vector y(m_Rotation.colY());
00328 Hep3Vector z(m_Rotation.colZ());
00329
00330 thetaX = kRad * x.theta();
00331 phiX = kRad * x.phi();
00332 thetaY = kRad * y.theta();
00333 phiY = kRad * y.phi();
00334 thetaZ = kRad * z.theta();
00335 phiZ = kRad * z.phi();
00336
00337
00338
00339
00340
00341
00342 return;
00343 }
00344
00345 int
00346 MucGeoGap::GuessStrip(const float x,
00347 const float y,
00348 const float ) const
00349 {
00350
00351 float x0, y0, z0;
00352 float xn, yn, zn;
00353 float dx;
00354 int n, iGuess;
00355
00356 n = m_StripNum - 1;
00357
00358 GetStrip(0)->GetCenterPos(x0, y0, z0);
00359 GetStrip(n)->GetCenterPos(xn, yn, zn);
00360
00361
00362
00363
00364
00365
00366 if(m_Orient == 0){
00367 dx = (yn-y0)/n;
00368 iGuess = int ((y-y0)/dx+0.5);
00369 return iGuess;
00370 }
00371 else{
00372 dx = (xn-x0)/n;
00373 iGuess = int ((x-x0)/dx+0.5);
00374 return iGuess;
00375 }
00376 }
00377
00378 HepPoint3D
00379 MucGeoGap::ProjectToGap(const HepPoint3D gPoint,
00380 const Hep3Vector gDirection) const
00381 {
00382
00383
00384 Hep3Vector gapVect = m_Rotation.colZ();
00385 HepPoint3D gCross(0,0,0);
00386
00387 HepPlane3D plane(gapVect, m_Center);
00388 MucGeometron geometron;
00389
00390 geometron.GetIntersectionLinePlane(gPoint, gDirection, plane, gCross);
00391
00392 return gCross;
00393 }
00394
00395 HepPoint3D
00396 MucGeoGap::ProjectToGapWithSigma(const HepPoint3D gPoint,
00397 const Hep3Vector gDirection,
00398 const HepPoint3D gPointSigma,
00399 const Hep3Vector gDirectionSigma,
00400 HepPoint3D &gCross,
00401 HepPoint3D &gCrossSigma) const
00402 {
00403 Hep3Vector gapVect = m_Rotation.colZ();
00404
00405 HepPlane3D plane(gapVect, m_Center);
00406 MucGeometron geometron;
00407 geometron.GetIntersectionLinePlaneWithSigma(gPoint, gDirection, gPointSigma, gDirectionSigma, plane, gCross, gCrossSigma);
00408
00409
00410 return gCross;
00411 }
00412
00413
00414 void
00415 MucGeoGap::ProjectToGapSurface(const HepPoint3D gPoint,
00416 const Hep3Vector gDirection,
00417 HepPoint3D& gCross1,
00418 HepPoint3D& gCross2 ) const
00419 {
00420
00421
00422 Hep3Vector gapVect = m_Rotation.colZ();
00423
00424 HepPlane3D planeInner(gapVect, m_SurfaceInner);
00425 HepPlane3D planeOuter(gapVect, m_SurfaceOuter);
00426 MucGeometron geometron;
00427
00428
00429 geometron.GetIntersectionLinePlane(gPoint, gDirection, planeInner, gCross1);
00430 geometron.GetIntersectionLinePlane(gPoint, gDirection, planeOuter, gCross2);
00431
00432 }
00433
00434 HepPoint3D
00435 MucGeoGap::ProjectToGapQuadLocal(const int part,
00436 const int orient,
00437 const float a,
00438 const float b,
00439 const float c,
00440 const int whichhalf,
00441 HepPoint3D& gCross1,
00442 HepPoint3D& gCross2) const
00443 {
00444 Hep3Vector gapVect, center_local;
00445 if(part == 1 && orient == 1) {gapVect = m_Rotation.colZ(); center_local = m_Center;}
00446 else if(part == 1 && orient == 0) {
00447 gapVect.setX(0); gapVect.setY(1); gapVect.setZ(0);
00448 center_local.setX(0); center_local.setY(m_Center.mag()); center_local.setZ(0);
00449 }
00450 else {
00451 gapVect.setX(1); gapVect.setY(0); gapVect.setZ(0);
00452 center_local.setX(m_Center.z()); center_local.setY(0); center_local.setZ(0);
00453 }
00454
00455
00456
00457
00458 HepPlane3D plane(gapVect, center_local);
00459
00460 MucGeometron geometron;
00461 geometron.GetIntersectionQuadPlaneLocal(part,orient, a, b, c, plane, gCross1, gCross2);
00462
00463
00464 }
00465
00466
00467 HepPoint3D
00468 MucGeoGap::ProjectToGap(const HepPoint3D gPoint,
00469 const float vy,
00470 const float y0,
00471 const float a,
00472 const float b,
00473 const float c,
00474 const int whichhalf,
00475 HepPoint3D& gCross1,
00476 HepPoint3D& gCross2) const
00477 {
00478
00479
00480 Hep3Vector gapVect = m_Rotation.colZ();
00481
00482 HepPlane3D plane(gapVect, m_Center);
00483 MucGeometron geometron;
00484
00485 geometron.GetIntersectionQuadPlane(gPoint, vy, y0, a, b, c, plane, gCross1, gCross2);
00486
00487 HepPoint3D localCross = TransformToGap(gCross1);
00488
00489 bool found = false;
00490 bool found1 = false;
00491 bool found2 = false;
00492 HepPoint3D temp;
00493 int good = 0;
00494
00495 if ( IsInGap(localCross.x(),
00496 localCross.y(),
00497 localCross.z()) ) {
00498 good = 1;
00499 found1 = true;
00500 }
00501
00502 localCross = TransformToGap(gCross2);
00503 if ( IsInGap(localCross.x(),
00504 localCross.y(),
00505 localCross.z()) ) {
00506 good = 2;
00507 found2 = true;
00508 }
00509 if(found1&&found2) {
00510 float center = b/(-2*a);
00511 if(whichhalf==2){
00512 if(gCross1.x()>center) good = 1;
00513 if(gCross2.x()>center) good = 2;
00514 }
00515 if(whichhalf==1){
00516 if(gCross1.x()<center) good = 1;
00517 if(gCross2.x()<center) good = 2;
00518 }
00519 }
00520 if(good == 2) {temp = gCross1; gCross1 = gCross2; gCross2 = temp;}
00521
00522
00523 return gCross1;
00524 }
00525
00526
00527 void
00528 MucGeoGap::ProjectToGapSurface(const HepPoint3D gPoint,
00529 const float vy,
00530 const float y0,
00531 const float a,
00532 const float b,
00533 const float c,
00534 const int whichhalf,
00535 HepPoint3D& gCross1,
00536 HepPoint3D& gCross2) const
00537 {
00538
00539
00540
00541 HepPoint3D cross1, cross2, cross3, cross4;
00542
00543 Hep3Vector gapVect = m_Rotation.colZ();
00544
00545 HepPlane3D planeInner(gapVect, m_SurfaceInner);
00546 HepPlane3D planeOuter(gapVect, m_SurfaceOuter);
00547
00548 MucGeometron geometron;
00549
00550 geometron.GetIntersectionQuadPlane(gPoint, vy, y0, a, b, c, planeInner, cross1, cross2);
00551 geometron.GetIntersectionQuadPlane(gPoint, vy, y0, a, b, c, planeOuter, cross3, cross4);
00552
00553 gCross1 = CompareIntersection(whichhalf, cross1, cross2, a, b, c);
00554 gCross2 = CompareIntersection(whichhalf, cross3, cross4, a, b, c);
00555
00556 }
00557
00558 HepPoint3D
00559 MucGeoGap::CompareIntersection(const int whichhalf, const HepPoint3D gCross1,
00560 const HepPoint3D gCross2, const float a,
00561 const float b, const float c )const
00562 {
00563 bool found = false;
00564 bool found1 = false;
00565 bool found2 = false;
00566 int good = 0;
00567
00568 HepPoint3D localCross = TransformToGap(gCross1);
00569 if ( IsInGap(localCross.x(),
00570 localCross.y(),
00571 0.0) ) {
00572 good = 1;
00573 found1 = true;
00574 }
00575
00576 localCross = TransformToGap(gCross2);
00577 if ( IsInGap(localCross.x(),
00578 localCross.y(),
00579 0.0) ) {
00580 good = 2;
00581 found2 = true;
00582 }
00583 if(found1&&found2) {
00584 float center = b/(-2*a);
00585 if(whichhalf==2){
00586 if(gCross1.x()>center) good = 1;
00587 if(gCross2.x()>center) good = 2;
00588 }
00589 if(whichhalf==1){
00590 if(gCross1.x()<center) good = 1;
00591 if(gCross2.x()<center) good = 2;
00592 }
00593 }
00594
00595 if(good == 1) return gCross1;
00596 else if(good == 2) return gCross2;
00597 else {
00598 HepPoint3D Empty(-9999,-9999,-9999);
00599 cout<<"in MucGeoGap:: both intersection position are bad!!!"<<endl; return Empty;
00600 }
00601
00602 }
00603
00604
00605 Hep3Vector
00606 MucGeoGap::RotateToGlobal(const Hep3Vector pVect) const
00607 {
00608
00609 return m_Rotation * pVect;
00610 }
00611
00612 Hep3Vector
00613 MucGeoGap:: RotateToGap(const Hep3Vector gVect) const
00614 {
00615
00616 return m_RotationT * gVect;
00617 }
00618
00619 HepPoint3D
00620 MucGeoGap::TransformToGlobal(const HepPoint3D pPoint) const
00621 {
00622
00623 return m_Rotation * pPoint + m_Translation;
00624 }
00625
00626 HepPoint3D
00627 MucGeoGap::TransformToGap(const HepPoint3D gPoint) const
00628 {
00629
00630 return m_RotationT * (gPoint - m_Translation);
00631 }
00632
00633 bool
00634 MucGeoGap::IsInGap(const float x,
00635 const float y,
00636 const float z) const
00637 {
00638
00639 return ( ( x > -0.5*m_XSize - kGapEdge ) && ( x < 0.5*m_XSize + kGapEdge ) &&
00640 ( y > -0.5*m_YSize - kGapEdge ) && ( y < 0.5*m_YSize + kGapEdge ) &&
00641 ( z > ( m_dzHighEdge - m_ZSize ) ) && ( z < m_dzHighEdge ) );
00642 }
00643
00644 MucGeoStrip*
00645 MucGeoGap::AddStrip(const int strip)
00646 {
00647
00648 MucGeoStrip *pStrip;
00649 MucGeoStrip *neighbor;
00650
00651 if( strip >= int(m_pMucGeoStrip.size()) ) {
00652 cout << " MucGeoGap::AddStrip strip number "
00653 << strip << " outside of expected range "
00654 << m_pMucGeoStrip.size()
00655 << endl;
00656 return 0;
00657 }
00658
00659 if( (strip + 1) > m_StripNum ) m_StripNum = strip + 1;
00660
00661 if(m_pMucGeoStrip[strip]) {
00662
00663 return m_pMucGeoStrip[strip];
00664 }
00665
00666 pStrip = new MucGeoStrip(strip, this);
00667
00668 m_pMucGeoStrip[strip] = pStrip;
00669 pStrip->SetLeftNeighbor(0L);
00670 pStrip->SetRightNeighbor(0L);
00671
00672
00673 if ( strip > 0 ) {
00674 neighbor = m_pMucGeoStrip[strip-1];
00675 if(neighbor) {
00676 neighbor->SetRightNeighbor(pStrip);
00677 pStrip->SetLeftNeighbor(neighbor);
00678 }
00679 }
00680
00681 return pStrip;
00682 }
00683
00684 ostream&
00685 operator << (ostream& s, const MucGeoGap& gap)
00686 {
00687
00688 float dx, dy, dz;
00689 HepPoint3D center;
00690 center = gap.GetCenter();
00691 gap.GetSize(dx,dy,dz);
00692
00693 s << " Identifier : " << gap.Part() << " " << gap.Seg() << " " << gap.Gap() << endl;
00694 s << " Strip number " << gap.GetStripNum() << endl;
00695 s << " Center position (" << center.x() << "," << center.y() << "," << center.z() << ")" << endl;
00696 s << " Size (" << dx << "," << dy << "," << dz << ")" << endl;
00697
00698 return s;
00699 }
00700
00701