StRoot: StBTofUtil/StBTofGeometry.cxx Source File

00001 
00002 /*******************************************************************
00003  *
00004  * $Id: StBTofGeometry.cxx,v 1.12 2011/07/27 16:15:12 geurts Exp $
00005  * 
00006  * Authors: Shuwei Ye, Xin Dong
00007  *******************************************************************
00008  *
00009  * Description: Collection of geometry classes for the TOF-MRPC
00010  *              initializes from GEANT geometry
00011  *
00012  *******************************************************************
00013  * $Log: StBTofGeometry.cxx,v $
00014  * Revision 1.12  2011/07/27 16:15:12  geurts
00015  * Alignment calibration modifications [Patrick Huck]:
00016  *  - added mAlignFile and SetAlignFile for use in StBTofMatchMaker
00017  *  - phi0, x0, z0 made mNTrays dependent
00018  *
00019  * Revision 1.11  2010/09/17 20:40:09  geurts
00020  * Protect Init() and InitFromStar() against non-initialized database/geant.
00021  * No immediate crash, but a LOG_ERROR instead.
00022  *
00023  * Revision 1.10  2010/08/09 18:45:36  geurts
00024  * Include methods in StBTofNode and StBTofGeometry that calculate local theta [Masa]
00025  *
00026  * Revision 1.9  2010/07/14 20:35:28  geurts
00027  * introduce switch to enable ideal MC geometry, without alignment updates. Default: disabled
00028  *
00029  * Revision 1.8  2010/05/25 22:09:44  geurts
00030  * improved database handling and reduced log output
00031  *
00032  * Revision 1.7  2010/04/03 02:00:53  dongx
00033  * X0 (radial offset) included in the tray alignment
00034  *
00035  * Revision 1.6  2009/09/15 00:17:27  dongx
00036  * Corrected the calculation for tray alignment parameters in X-Y
00037  *
00038  * Revision 1.5  2009/08/25 15:41:29  fine
00039  * fix the compilation issues under SL5_64_bits  gcc 4.3.2
00040  *
00041  * Revision 1.4  2009/03/18 14:18:18  dongx
00042  * - Optimized the geometry initialization function, reduced the CPU time use
00043  * - Optimized the HelixCrossCellIds() function, now doing the tray fast projection to reduce the loop
00044  *
00045  * Revision 1.3  2009/02/13 00:00:56  dongx
00046  * Tray geometry alignment implemented.
00047  *
00048  * Revision 1.2  2009/02/12 01:45:57  dongx
00049  * Clean up
00050  *
00051  * Revision 1.1  2009/02/02 21:56:54  dongx
00052  * first release - Barrel geometry
00053  *
00054  *******************************************************************/
00055 #include "Stiostream.h"
00056 #include <math.h>
00057 #include <vector>
00058 #include <string>
00059 #include <stdlib.h>
00060 #include <stdio.h>
00061 #include "tables/St_tofGeomAlign_Table.h"
00062 
00063 #include "StBTofGeometry.h"
00064 #include "TFile.h"
00065 #include "TDataSet.h"
00066 #include "TDataSetIter.h"
00067 #include "StMaker.h"
00068 //#include "TMemStat.h"
00069 #include "StMessMgr.h"
00070 
00071 //#include "Debugger.h"
00072 
00074 //
00075 // group of classes for BTof Geometry:
00076 //
00077 //    StBTofGeometry, StBTofNode
00078 //    StBTofGeomTray, StBTofGeomSensor
00079 //
00080 // Usage:
00081 //
00082 //   StBTofGeometry* geo = new StBTofGeometry("btof","btof geometry");
00083 //    geo->Init(TVolume *starHall);
00084 //
00085 // ---------------------------------------------------------------------------
00086 //
00087 // StBTofNode
00088 // ==============
00089 //
00091 
00092 #ifdef __ROOT__      
00093 ClassImp(StBTofNode)
00094 #endif
00095 
00096 Bool_t StBTofNode::mDebug = kFALSE;
00097 Double_t const StBTofGeomSensor::mSensorDy = 10.35;   // Actual module length;
00098 const char* StBTofGeometry::sectorPref = "BSEC";
00099 const char* StBTofGeometry::trayPref   = "BTRA";
00100 const char* StBTofGeometry::senPref    = "BRMD";
00101 
00102 //_____________________________________________________________________________
00103 StBTofNode::StBTofNode(TVolumeView *element, TVolumeView *top, StThreeVectorD *align, TVolumePosition *pos)
00104   : fView(element), pView(new TVolumePosition(*pos)), mMasterNode(top), mTransFlag(kFALSE)
00105 {
00106    if(align) {
00107      mAlign[0] = align->x();  mAlign[1] = align->y();  mAlign[2] = align->z();
00108    } else {
00109      mAlign[0] = 0.0;   mAlign[1] = 0.0;   mAlign[2] = 0.0;
00110    }
00111 
00112    UpdateMatrix();
00113    BuildMembers();
00114 }
00115 
00116 //_____________________________________________________________________________
00117 StBTofNode::~StBTofNode()
00118 { 
00119   fView = 0;
00120   pView = 0;
00121   mMasterNode = 0;
00122   mTransFlag = kFALSE;
00123 }
00124 
00125 //_____________________________________________________________________________
00126 void  StBTofNode::UpdateMatrix()
00127 {
00128    //
00129    //Update the Translation and RotateMatrix between local and master
00130    //  and store them as members: mTransMRS, mRotMRS
00131    //while TNode stores them locally to share among all TNode objects,
00132    //  thus they may be changed afterward !
00133    //
00134 
00135    mTransFlag = kFALSE;
00136    CalcMatrix(this, mAlign, mTransMRS, mRotMRS);
00137    mTransFlag = kTRUE;
00138 
00139 }
00140 
00141 //_____________________________________________________________________________
00142 void StBTofNode::CalcMatrix(StBTofNode* son, Double_t* align,
00143               Double_t* trans, Double_t* rot, StBTofNode* mother)
00144 {
00145    //
00146    //Translation vector and Rotation matrix from TNode "mother" to "son"
00147    //
00148    // 1. Transformation: mother(x,y,z) ==> son(x',y',z')
00149    //
00150    //    [x']    [ m[0]  m[1]  m[2] ]   [ [x]   [dx] ]
00151    //    [y']  = [ m[3]  m[4]  m[5] ] . [ [y] - [dy] ]
00152    //    [z']    [ m[6]  m[7]  m[8] ]   [ [z]   [dz] ]
00153    //
00154    //       Transpose([x',y',z'])       = M . Transpose([x-dx,y-dy,z-dz])
00155    //  and
00156    //       Transpose([x-dx,y-dy,z-dz]) = inverse(M) . Transpose([x',y',z'])
00157    //
00158    //     where transpose(M) = inverse(M)
00159    //
00160    // 2. In the representation of vector:
00161    //
00162    // Vectors: V' = x' I' + y' J' + z' K'
00163    //          V  = x  I  + y  J  + z  K   =  V' + dV
00164    //         dV  = dx I  + dy J  + dz K
00165    //
00166    //    V' . I' = x' = (V-dV) . I'
00167    //            = (x-dx) (I.I') + (y-dy) (J.I') + (z-dz) (K.I')
00168    //            = (x-dx) m[0]   + (y-dy) m[1]   + (z-dz) m[2]
00169    //
00170    //   so, the projections of unit vector I' on I, J, K
00171    //       yield m[0], m[1], m[2]
00172    //
00173    //********************************************************************
00174 
00175 
00176    Double_t xl[3], xm[3];
00177 
00178    //  son->TNode::UpdateMatrix();
00179    //   son->GetpView()->UpdateMatrix();
00180    //   LOG_INFO << "StBTofNode::CalcMatrix" << endm;
00181 
00182    xl[0] = 0;  xl[1] = 0;  xl[2] = 0;
00183    ConvertPos(son,xl, mother,xm);
00184    trans[0] = xm[0];
00185    trans[1] = xm[1];
00186    trans[2] = xm[2];
00187 
00188    xl[0] = 1;  xl[1] = 0;  xl[2] = 0;
00189    ConvertPos(son,xl, mother,xm);
00190    rot[0] = xm[0]-trans[0];
00191    rot[1] = xm[1]-trans[1];
00192    rot[2] = xm[2]-trans[2];
00193 
00194    xl[0] = 0;  xl[1] = 1;  xl[2] = 0;
00195    ConvertPos(son,xl, mother,xm);
00196    rot[3] = xm[0]-trans[0];
00197    rot[4] = xm[1]-trans[1];
00198    rot[5] = xm[2]-trans[2];
00199 
00200    xl[0] = 0;  xl[1] = 0;  xl[2] = 1;
00201    ConvertPos(son,xl, mother,xm);
00202    rot[6] = xm[0]-trans[0];
00203    rot[7] = xm[1]-trans[1];
00204    rot[8] = xm[2]-trans[2];
00205 
00206    // Alignment parameters;
00207    trans[0] += align[0];
00208    trans[1] += align[1];
00209    trans[2] += align[2];   
00210 
00211    return;
00212 }
00213 
00214 //_____________________________________________________________________________
00215 void StBTofNode::ConvertPos(
00216                      StBTofNode* from, const Double_t* pos_from,
00217             StBTofNode* to,         Double_t* pos_to)
00218 {
00219    if (to==0) from->Local2Master(pos_from,pos_to);
00220    else {
00221       Double_t xg[3];
00222       from->Local2Master(pos_from,xg);
00223       to->Master2Local(xg,pos_to);
00224    }
00225 
00226    return;
00227 }
00228 
00229 //_____________________________________________________________________________
00230 void  StBTofNode::BuildMembers()
00231 {
00232    //
00233    //Build date members: mCenter{Rxy,Eta,Phi}, m{Eta,Phi}{Min,Max}
00234    //
00235 
00236    Double_t xl[3];
00237    Double_t xg[3];
00238    TBRIK *brik = dynamic_cast<TBRIK*>(GetShape());
00239    //   LOG_INFO << " Get shape ready" << endm;
00240    if(!brik) { LOG_INFO << " no brik " << endm; return; }
00241    //   Double_t dx = brik->GetDx();
00242    Double_t dy = brik->GetDy();
00243    Double_t dz = brik->GetDz();
00244 
00245    //   LOG_INFO << " size = " <<dx << " " << dy << " " << dz << endm;
00246 /*
00247    //center point
00248    Double_t xc[3];
00249    xc[0] = 0;  xc[1] = 0;  xc[2] = 0;
00250    Local2Master(xc, xg);
00251    StThreeVectorD center(xg[0], xg[1], xg[2]);
00252    mCenterRxy = center.perp();
00253    mCenterEta = center.pseudoRapidity();
00254    mCenterPhi = center.phi();
00255 */
00256 
00257    // -dz
00258    xl[0] = 0;  xl[1] = 0;  xl[2] = -dz;
00259    Local2Master(xl, xg);
00260    mEtaMin = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
00261 
00262    // +dz
00263    xl[0] = 0;  xl[1] = 0;  xl[2] = dz;
00264    Local2Master(xl, xg);
00265    mEtaMax = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
00266 
00267    // -dy
00268    xl[0] = 0;  xl[1] = -dy;  xl[2] = 0;
00269    Local2Master(xl, xg);
00270    mPhiMin = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
00271 
00272    // +dy
00273    xl[0] = 0;  xl[1] =  dy;  xl[2] = 0;
00274    Local2Master(xl, xg);
00275    mPhiMax = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
00276 
00277    if (mEtaMin>mEtaMax) {
00278       Double_t tmp = mEtaMin;
00279       mEtaMin = mEtaMax;
00280       mEtaMax = tmp;
00281    }
00282 
00283    if (mPhiMin>mPhiMax) {
00284       Double_t tmp = mPhiMin;
00285       mPhiMin = mPhiMax;
00286       mPhiMax = tmp;
00287    }
00288 
00289    return;
00290 }
00291 
00292 //_____________________________________________________________________________
00293 void StBTofNode::Local2Master(const Double_t* local, Double_t* master)
00294 {
00295    //
00296    //Transform local coordinate into global coordinate
00297    //
00298    //  pView->UpdateMatrix();
00299    if (!mTransFlag) {
00300      if (!mMasterNode) { LOG_INFO << " no Master! " << endm; return;}
00301      if (pView) {
00302        pView->Local2Master(local, master);
00303      } else {
00304        TVolumeView *son = GetfView();
00305        TVolumeView *mrs = GetTopNode();
00306      
00307        TVolumePosition *pos = 0;
00308        pos = son->Local2Master(son, mrs);
00309        pos->Local2Master(local, master);
00310        delete pos;
00311      }
00312      return;
00313    }
00314 
00315    //mTransFlag==kTRUE, i.e. StBTofGeomNode::UpdateMatrix() invoked already
00316    Double_t x, y, z;
00317    x = local[0];
00318    y = local[1];
00319    z = local[2];
00320 
00321    master[0] = mTransMRS[0] + mRotMRS[0]*x + mRotMRS[3]*y + mRotMRS[6]*z;
00322    master[1] = mTransMRS[1] + mRotMRS[1]*x + mRotMRS[4]*y + mRotMRS[7]*z;
00323    master[2] = mTransMRS[2] + mRotMRS[2]*x + mRotMRS[5]*y + mRotMRS[8]*z;
00324 
00325 }
00326 
00327 //_____________________________________________________________________________
00328 void StBTofNode::Master2Local(const Double_t* master, Double_t* local)
00329 {
00330    //
00331    //Transform global coordinate into local coordinate
00332    //
00333    //  pView->UpdateMatrix();
00334    //   pView->Master2Local(master, local);
00335    if (!mTransFlag) {
00336      LOG_INFO << " and No TVolumePosition::Master2Local is wrong, so do nothing" << endm;
00337       return;
00338    }
00339 
00340    //mTransFlag==kTRUE, i.e. StBTofGeomNode::UpdateMatrix() invoked already
00341    Double_t x, y, z;
00342    x = master[0] - mTransMRS[0];
00343    y = master[1] - mTransMRS[1];
00344    z = master[2] - mTransMRS[2];
00345 
00346    local[0] = mRotMRS[0]*x + mRotMRS[1]*y + mRotMRS[2]*z;
00347    local[1] = mRotMRS[3]*x + mRotMRS[4]*y + mRotMRS[5]*z;
00348    local[2] = mRotMRS[6]*x + mRotMRS[7]*y + mRotMRS[8]*z;
00349 }
00350 
00351 //_____________________________________________________________________________
00352 StThreeVectorD StBTofNode::YZPlaneNormal()
00353 {
00354    //
00355    //Calculate the vector unit of normal to YZ-plane
00356    // i.e. the global representation of local unit vector (1,0,0)
00357    //
00358 
00359    Double_t ux[3], nx[3];
00360 
00361    ux[0] = 1;  ux[1] = 0;  ux[2] = 0;
00362    Local2Master(ux,nx);
00363 
00364    nx[0] -= mTransMRS[0];
00365    nx[1] -= mTransMRS[1];
00366    nx[2] -= mTransMRS[2];
00367 
00368    return StThreeVectorD(nx[0],nx[1],nx[2]);
00369 }
00370 
00371 //_____________________________________________________________________________
00372 StThreeVectorD  StBTofNode::GetCenterPosition()
00373 const
00374 {
00375    //
00376    //Return the global representation of this node(TBRIk-shape) center
00377    //
00378 
00379    Double_t xg[3];
00380    xg[0] = mTransMRS[0];
00381    xg[1] = mTransMRS[1];
00382    xg[2] = mTransMRS[2];
00383 
00384    return StThreeVectorD(xg[0],xg[1],xg[2]);
00385 }
00386 
00387 //_____________________________________________________________________________
00388 Bool_t StBTofNode::IsLocalPointIn(const Double_t x, const Double_t y, 
00389                                       const Double_t z)
00390 {
00391    TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
00392    Double_t dx = brik->GetDx();
00393    Double_t dy = brik->GetDy();
00394    Double_t dz = brik->GetDz();
00395    Bool_t ret = -dx<x && x<dx && -dy<y && y<dy && -dz<z && z<dz;
00396 
00397    return ret;
00398 }
00399 
00400 //_____________________________________________________________________________
00401 Bool_t  StBTofNode::IsGlobalPointIn(const StThreeVectorD &global)
00402 {
00403    Double_t xl[3], xg[3];
00404    xg[0] = global.x();
00405    xg[1] = global.y();
00406    xg[2] = global.z();
00407    Master2Local(xg, xl);
00408    Bool_t ret = IsLocalPointIn(xl[0], xl[1], xl[2]);
00409 
00410    return ret;
00411 }
00412 
00413 //_____________________________________________________________________________
00414 Bool_t StBTofNode::HelixCross(const StHelixD &helix, Double_t &pathLen,
00415                                   StThreeVectorD &cross)
00416 {
00417    //
00418    // check if helix go through this node(TBRIK)
00419    //  and return the path length of helix before crossing this node
00420    //
00421    //   static const Float_t MaxPathLength = 1000.; //Maximum path length
00422    Float_t MaxPathLength = 1000.;
00423 
00424    Bool_t ret = kFALSE;
00425    pathLen = 0;
00426 
00427    //
00428    // Get the normal to the YZ-plane
00429    //
00430    StThreeVectorD planeNormal = YZPlaneNormal();
00431 
00432    //
00433    // Get the center position
00434    //
00435    StThreeVectorD centerPos = GetCenterPosition();
00436 
00437    //
00438    // Find the intersection point between the helix and the cell plane
00439    //
00440    pathLen = helix.pathLength(centerPos,planeNormal);
00441    if ( pathLen>0 && pathLen<MaxPathLength ) {
00442       cross = helix.at(pathLen);
00443       //
00444       // Check if the intersected point is really in the cell
00445       //
00446       ret = IsGlobalPointIn(cross);
00447    }
00448    return ret;
00449 }
00450 
00451 //_____________________________________________________________________________
00452 Bool_t StBTofNode::HelixCross(const StHelixD &helix, Double_t &pathLen,
00453                                   StThreeVectorD &cross, Double_t &theta)
00454 {
00455    //
00456    // check if helix go through this node(TBRIK)
00457    //  and return the path length of helix before crossing this node
00458    //
00459    //   static const Float_t MaxPathLength = 1000.; //Maximum path length
00460    Float_t MaxPathLength = 1000.;
00461 
00462    Bool_t ret = kFALSE;
00463    pathLen = 0;
00464 
00465    //
00466    // Get the normal to the YZ-plane
00467    //
00468    StThreeVectorD planeNormal = YZPlaneNormal();
00469 
00470    //
00471    // Get the center position
00472    //
00473    StThreeVectorD centerPos = GetCenterPosition();
00474 
00475    //
00476    // Find the intersection point between the helix and the cell plane
00477    //
00478    pathLen = helix.pathLength(centerPos,planeNormal);
00479    if ( pathLen>0 && pathLen<MaxPathLength ) {
00480            cross = helix.at(pathLen);
00481            theta = planeNormal.angle(helix.cat(pathLen));
00482       //
00483       // Check if the intersected point is really in the cell
00484       //
00485       ret = IsGlobalPointIn(cross);
00486    }
00487    return ret;
00488 }
00489 
00490 //_____________________________________________________________________________
00491 void StBTofNode::Print(Option_t *opt) const
00492 {
00493    TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
00494    LOG_INFO << "Name=" << GetName() << "\t TBRIK-dimension=" << brik->GetDx()
00495         << " : " << brik->GetDy() << " : " << brik->GetDz()
00496         // << "\n Center Rxy:Eta:Phi=" << mCenterRxy << " : "
00497         // << mCenterEta << " : " << mCenterPhi
00498         << "\n EtaRange=" << mEtaMin << " : " << mEtaMax
00499         << "\t PhiRange=" << mPhiMin << " : " << mPhiMax
00500         << endm;
00501    LOG_INFO <<"trans[0-2]=" << mTransMRS[0] <<" " <<mTransMRS[1]
00502                                          <<" " <<mTransMRS[2]
00503         <<"\nmRotMRS[0-2]=" <<mRotMRS[0] <<" " <<mRotMRS[1] <<" " <<mRotMRS[2]
00504         <<"\nmRotMRS[3-5]=" <<mRotMRS[3] <<" " <<mRotMRS[4] <<" " <<mRotMRS[5]
00505         <<"\nmRotMRS[6-8]=" <<mRotMRS[6] <<" " <<mRotMRS[7] <<" " <<mRotMRS[8]
00506         <<endm;
00507 }
00508 
00510 //
00511 // StBTofGeomTray
00512 // ==============
00513 //
00515 
00516 #ifdef __ROOT__      
00517 ClassImp(StBTofGeomTray)
00518 #endif
00519 
00520 Bool_t StBTofGeomTray::mDebug = kFALSE;
00521 
00522 //_____________________________________________________________________________
00523 StBTofGeomTray::StBTofGeomTray(const Int_t ibtoh, TVolumeView *sector, TVolumeView *top, StThreeVectorD *align, TVolumePosition *pos)
00524   : StBTofNode((TVolumeView *)sector->First(), top, align, pos)
00525 {
00526   mSectorsInBTOH = top->GetListSize()/2;
00527   mBTOHIndex = ibtoh + 1;
00528   mTrayIndex = ibtoh * mSectorsInBTOH + sector->GetPosition()->GetId();
00529 }
00530 
00531 //_____________________________________________________________________________
00532 StBTofGeomTray::~StBTofGeomTray()
00533 {
00534   mBTOHIndex = 0;
00535   mTrayIndex = 0;
00536 }
00537 //_____________________________________________________________________________
00538 void StBTofGeomTray::Print(const Option_t *opt) const
00539 {
00540    LOG_INFO << "StBTofGeomTray, tray#=" << mTrayIndex << endm;
00541    StBTofNode::Print(opt);
00542 }
00543 
00544 
00546 //
00547 // StBTofGeomSensor
00548 // ================
00549 //
00551 
00552 #ifdef __ROOT__      
00553 ClassImp(StBTofGeomSensor)
00554 #endif
00555 
00556 Bool_t StBTofGeomSensor::mDebug = kFALSE;
00557 
00558 //_____________________________________________________________________________
00559 StBTofGeomSensor::StBTofGeomSensor(TVolumeView *element, TVolumeView *top, StThreeVectorD *align, TVolumePosition *pos) 
00560   : StBTofNode(element, top, align, pos)
00561 {
00562    mModuleIndex = element->GetPosition()->GetId();
00563    CreateGeomCells();
00564 }
00565 
00566 //_____________________________________________________________________________
00567 StBTofGeomSensor::~StBTofGeomSensor()
00568 {
00569   mModuleIndex = 0;
00570 }
00571 //_____________________________________________________________________________
00572 void StBTofGeomSensor::CreateGeomCells()
00573 {
00574    //
00575    //Divide this sensor to creat cells
00576    //
00577 
00578    //
00579    // change the size according to the real cells
00580    //    mSensorDy -- defined by myself
00581    //
00582    Double_t sensor_dy = mSensorDy;
00583    Double_t cell_width = 2*sensor_dy/mCells;
00584 
00585    for (int i=0; i<=mCells; i++) mCellY[i] = cell_width*i - sensor_dy;
00586 
00587 }
00588 
00589 //_____________________________________________________________________________
00590 Double_t StBTofGeomSensor::GetCellYMin(const Int_t icell)
00591 const
00592 {
00593    Double_t ret = 0;
00594    if (icell<=0 || icell>mCells) {
00595       LOG_INFO << "cell#=" << icell <<" is out range=[0," << mCells << "]"
00596            << endm;
00597    } else ret = mCellY[icell-1];
00598 
00599    return ret;
00600 }
00601 
00602 //_____________________________________________________________________________
00603 Double_t StBTofGeomSensor::GetCellYMax(const Int_t icell)
00604 const
00605 {
00606    Double_t ret = 0;
00607    if (icell<=0 || icell>mCells) {
00608       LOG_INFO << "cell#=" << icell <<" is out range=[0," << mCells << "]"
00609            << endm;
00610    } else ret = mCellY[icell];
00611 
00612    return ret;
00613 }
00614 
00615 //_____________________________________________________________________________
00616 StThreeVectorD StBTofGeomSensor::GetCellPosition(const Int_t icell)
00617 {
00618    //
00619    // Get the center position of cell in this sensor
00620    //
00621 
00622    static const char* where = "StBTofGeomSensor::GetCellPosition";
00623    /*
00624    TBRIK *thisBrik = dynamic_cast<TBRIK*> (GetShape());
00625    Double_t sensor_dy  = thisBrik->GetDy();
00626    */
00627 
00628    // change the size according to the real cells
00629    Double_t sensor_dy = mSensorDy;
00630    Double_t cell_dy    = 2*sensor_dy/mCells;
00631 
00632    Double_t xl[3], xg[3];
00633    if (icell>=1 && icell<=mCells) {
00634       xl[0] = 0;
00635       xl[1] = (icell*2-1)*cell_dy - sensor_dy;
00636       xl[2] = 0;
00637       Local2Master(xl,xg);
00638    } else {  //invalid cell, return (0.,0.,0.)
00639       LOG_INFO << "Warning in " << where <<" Invalid cell# =" << icell << endm;
00640       xg[0] = 0.;
00641       xg[1] = 0.;
00642       xg[2] = 0.;
00643    }
00644 
00645    return StThreeVectorD(xg[0],xg[1],xg[2]);  
00646 }
00647 
00648 
00649 //_____________________________________________________________________________
00650 Int_t StBTofGeomSensor::FindCellIndex(const Double_t* local)
00651 {
00652    //
00653    //Look up the cell the local point in
00654    //
00655 
00656    Int_t icell=-1;
00657    if ( IsLocalPointIn(local[0],local[1],local[2]) ) {
00658     
00659       for (Int_t i=0; i<mCells; i++) {
00660          if (mCellY[i]<= local[1] && local[1]<=mCellY[i+1]) {
00661             icell = i+1;
00662             break;
00663             
00664          }
00665       }
00666    }
00667 
00668    return icell;
00669 }
00670 
00671 //_____________________________________________________________________________
00672 void StBTofGeomSensor::Print(const Option_t *opt) const
00673 {
00674    LOG_INFO << "StBTofGeomSensor, module#=" << mModuleIndex << endm;
00675    StBTofNode::Print(opt);
00676 
00677    LOG_INFO << " Cells=" << mCells << "\t Y range for cells=\n";
00678    for (Int_t i=0; i<=mCells; i++) LOG_INFO << " : " << mCellY[i];
00679    LOG_INFO << endm;
00680 }
00681 
00682 
00684 //
00685 // StBTofGeometry
00686 // ==============
00687 //
00689 
00690 StBTofGeometry *gBTofGeometry = 0;
00691 static const Int_t CELLSINMODULE = 6;
00692 
00693 #ifdef __ROOT__      
00694 ClassImp(StBTofGeometry)
00695 #endif
00696 
00697 Bool_t StBTofGeometry::mDebug = kFALSE;
00698 
00699 //_____________________________________________________________________________
00700 StBTofGeometry::StBTofGeometry(const char* name, const char* title)
00701   : TNamed(name,title)
00702 {
00703    mCellsInModule  = StBTofGeomSensor::GetCells();
00704    mModulesInTray  = 0;
00705    mNValidTrays    = 0;
00706    mRootFile       = 0;
00707    mInitFlag       = kFALSE;
00708    mTopNode        = 0;
00709    mStarHall       = 0;
00710    mIsMC           = kFALSE;
00711    SetAlignFile("");
00712 
00713    for(int i=0;i<mNTrays;i++) {
00714      mBTofTray[i] = 0;
00715      for(int j=0;j<mNModules;j++) {
00716        mBTofSensor[i][j] = 0;
00717      }
00718    }
00719 
00720    //
00721    //We only need one instance of StBTofGeometry
00722    //
00723    if (gBTofGeometry) {
00724       LOG_INFO << "Warning !! There is already StBTofGeometry at pointer="
00725            << (void*)gBTofGeometry << ", so it is deleted"
00726            << endm;
00727       delete gBTofGeometry;
00728    }
00729    gBTofGeometry = this;
00730 }
00731 
00732 //_____________________________________________________________________________
00733 StBTofGeometry::~StBTofGeometry()
00734 {
00735    LOG_INFO << "Warning !! StBTofGeometry at pointer =" << (void*)gBTofGeometry
00736         << " is deleted" << endm;
00737    gBTofGeometry = 0;
00738 
00739    for(int i=0;i<mNTrays;i++) {
00740      if(mBTofTray[i]) delete mBTofTray[i];
00741      mBTofTray[i] = 0;
00742      for(int j=0;j<mNModules;j++) {
00743        if(mBTofSensor[i][j]) delete mBTofSensor[i][j];
00744        mBTofSensor[i][j] = 0;
00745      }
00746    }
00747    
00748 }
00749 
00750 //_____________________________________________________________________________
00751 void StBTofGeometry::Init(StMaker *maker, TVolume *starHall)
00752 {
00753    //
00754    //Define geometry parameters and establish the geometry
00755    //  
00756    if(maker->Debug()) DebugOn();
00757 
00758    // retrieve align parameters  -- need to use db
00759    for(int i=0;i<mNTrays;i++) {
00760      mTrayX0[i] = 0.0;
00761      mTrayY0[i] = 0.0;
00762      mTrayZ0[i] = 0.0;
00763    }
00764 
00765    double phi0[mNTrays], x0[mNTrays], z0[mNTrays];
00766    for(int i=0;i<mNTrays;i++) {
00767      phi0[i] = 0.0, x0[i] = 0.0; z0[i] = 0.0;
00768    }
00769 
00770    // If not MC input, load the alignment parameters from the database; otherwise ignore.
00771    if (mIsMC) {
00772      LOG_INFO << "[StBTofGeometry] detected MC-mode: ignore alignment corrections" << endm;
00773    } else if (strcmp(mAlignFile.c_str(),"")!=0) {
00774      LOG_INFO << "[StBTofGeometry] receiving alignment parameters from input files" << endm;
00775      ifstream inData;
00776      inData.open(mAlignFile.c_str());
00777      if(inData.good()) {
00778        for(int i=0;i<mNTrays;i++) {
00779          inData >> phi0[i] >> z0[i] >> x0[i];
00780        }
00781      } else {
00782        LOG_WARN << " Bad input file ! Use ideal geometry! " << endm;
00783      }
00784      inData.close();
00785    } else {
00786      LOG_INFO << "[StBTofGeometry] retrieving geometry alignment parameters from database" << endm;
00787      TDataSet *mDbTOFDataSet = maker->GetDataBase("Calibrations/tof/tofGeomAlign");
00788      if (!mDbTOFDataSet) {
00789        LOG_WARN << "[StBTofGeometry] unable to find Calibrations/tof/tofGeomAlign! Use ideal geometry!" << endm;
00790      } else {
00791        St_tofGeomAlign* tofGeomAlign = static_cast<St_tofGeomAlign*>(mDbTOFDataSet->Find("tofGeomAlign"));
00792        if(!tofGeomAlign) {
00793          LOG_WARN << "Unable to get tof geometry align parameter! Use ideal geometry!" << endm;
00794        } else {
00795          tofGeomAlign_st* geomAlign = static_cast<tofGeomAlign_st*>(tofGeomAlign->GetArray());
00796      
00797          for (Int_t i=0;i<mNTrays;i++) {
00798            phi0[i] = geomAlign[i].phi0;
00799            x0[i]   = geomAlign[i].x0;
00800            z0[i]   = geomAlign[i].z0;
00801          }
00802        }
00803      }
00804    }
00805 
00806    for(int i=0;i<mNTrays;i++) {
00807      double phi;
00808      if(i<60) {
00809        phi = 72 - i*6;   // phi angle of tray Id = i+1, west
00810        double cs = TMath::Cos(phi*TMath::Pi()/180.);
00811        double ss = TMath::Sin(phi*TMath::Pi()/180.);
00812        mTrayX0[i] = phi0[i]*ss + x0[i]*cs;
00813        mTrayY0[i] = -phi0[i]*cs + x0[i]*ss;
00814      } else {
00815        phi = 108 + (i-60)*6;   // phi angle of tray Id = i+1, east
00816        double cs = TMath::Cos(phi*TMath::Pi()/180.);
00817        double ss = TMath::Sin(phi*TMath::Pi()/180.);
00818        mTrayX0[i] = -phi0[i]*ss + x0[i]*cs;
00819        mTrayY0[i] = phi0[i]*cs + x0[i]*ss;
00820      }
00821      mTrayZ0[i] = z0[i];
00822 
00823      if(maker->Debug()) {
00824        LOG_DEBUG << " Tray # = " << i+1 << " Align parameters " << mTrayX0[i] << " " << mTrayY0[i] << " " << mTrayZ0[i] << endm;
00825      }
00826    }
00827 
00828    InitFromStar(starHall);
00829    mStarHall = starHall;
00830 }
00831 //_____________________________________________________________________________
00832 void StBTofGeometry::InitFromStar(TVolume *starHall)
00833 {
00834   // Initialize TOFr geometry from STAR geometry
00835   //     BTofConf   --     0     tray_BTof   (default)
00836   //                       1     full_BTof
00837 
00838   if (!starHall) {
00839     LOG_ERROR << "[StBTofGeometry] No STAR Hall volume defined" << endm;
00840     return;
00841   }
00842 
00843   // Loop over the STAR geometry and mark the volume needed
00844   TDataSetIter volume(starHall,0);
00845 
00846   TVolume *starDetectorElement = 0;
00847   while ( (starDetectorElement = ( TVolume *)volume()) )
00848     {
00849       //      const char *elementName = starDetectorElement->GetName();
00850       //      Bool_t found = ! ( strcmp(elementName,tofElements[0]) && strcmp(elementName,tofElements[1]) && strcmp(elementName,tofElements[2]) );
00851       Bool_t found = ( IsBSEC(starDetectorElement) || IsBTRA(starDetectorElement) || IsBRMD(starDetectorElement) );
00852       if (found) {
00853 
00854         starDetectorElement->SetVisibility(TVolume::kBothVisible);
00855         starDetectorElement->Mark();
00856         if (starDetectorElement->GetLineColor()==1 || starDetectorElement->GetLineColor()==7) 
00857           starDetectorElement->SetLineColor(14);
00858         
00859       } else {
00860         
00861         starDetectorElement->UnMark();
00862         starDetectorElement->SetVisibility(TVolume::kThisUnvisible);
00863 
00864       }
00865     }
00866 
00867   starHall->SetVisibility(TVolume::kBothVisible);
00868   mTopNode = new TVolumeView(*starHall,10); 
00869 
00870   mSectorsInBTOH = mTopNode->GetListSize()/2;    // # of sectors in one half
00871 
00873   // save the sensors and trays
00875   if(!mTopNode) {
00876     LOG_WARN << " No Top Node for Tof Geometry! " << endm;
00877     return;
00878   }
00879   LOG_INFO << "[StBTofGeometry] # of sectors = " << mTopNode->GetListSize() << endm;
00880 
00881   TVolumePosition *transPos = 0;
00882 
00883 //  TDataSetIter nextDet(mTopNode, 0);
00884   TVolumeViewIter nextDet((TVolumeView *)mTopNode, 0);
00885   TVolumeView *secVolume = 0;
00886   TVolumeView *detVolume = 0;
00887   Int_t ibtoh = 0;
00888   Int_t isec = 0;
00889   Int_t isensor = 0;
00890   mNValidTrays = 0;
00891   mModulesInTray = 0;
00892   StThreeVectorD *align = 0;
00893   while ( (detVolume = (TVolumeView *)nextDet()) ) {
00894 
00895     if(strcmp(detVolume->GetName(), sectorPref)==0) {  // sector volume
00896       isec++;
00897       secVolume = (TVolumeView *)detVolume;
00898       continue;   // sector continue;
00899     }
00900     if(strcmp(detVolume->GetName(), trayPref)==0) {   // tray volume
00901       if ( isec>60 ) ibtoh = 1;
00902       int sectorsInBTOH = mTopNode->GetListSize()/2;
00903       int trayIndex = ibtoh * sectorsInBTOH + secVolume->GetPosition()->GetId(); // secVolume
00904       LOG_DEBUG << " Tray # " << trayIndex << " has # of modules = " << detVolume->GetListSize() << endm;
00905       isensor = 0;   // clear for this tray
00906       if(detVolume->GetListSize()) {   // valid tray
00907 
00908         int itray = trayIndex - 1;
00909 
00910         if(align) delete align;    align = 0;
00911         if(trayIndex<=60) align = new StThreeVectorD(mTrayX0[itray], mTrayY0[itray], mTrayZ0[itray]);
00912         else align = new StThreeVectorD(mTrayX0[itray], mTrayY0[itray], -mTrayZ0[itray]);  // thus z0 will be the distance between the tray end to the TPC central membrane
00913 
00914         transPos = nextDet[0];
00915 
00916         mNValidTrays++;
00917 
00918         mBTofTray[mNValidTrays-1] = new StBTofGeomTray(ibtoh, secVolume, mTopNode, align, transPos);
00919         delete transPos;  transPos = 0;
00920 
00921         if(mDebug) {
00922           LOG_DEBUG << "   Initialize and save tray # " << mBTofTray[mNValidTrays-1]->Index() << " with " << detVolume->GetListSize() << " modules" << endm;
00923           LOG_DEBUG << "   alignment parameters \t" << mTrayX0[itray] << " " <<  mTrayY0[itray] << " " <<  mTrayZ0[itray] << endm;
00924           mBTofTray[mNValidTrays-1]->Print();
00925         }
00926 
00927       }
00928     }
00929     if(strcmp(detVolume->GetName(), senPref)==0) {   // module volume
00930 
00931       transPos = nextDet[0];
00932 
00933       mBTofSensor[mNValidTrays-1][isensor] = new StBTofGeomSensor(detVolume, mTopNode, align, transPos);
00934       delete transPos;   transPos = 0;
00935 
00936       if(mDebug) mBTofSensor[mNValidTrays-1][isensor]->Print();
00937 
00938       isensor++;
00939       if(isensor>mModulesInTray) mModulesInTray = isensor;
00940     }
00941   }
00942 
00943   mBTofConf = 0;
00944   if (mNValidTrays==120) mBTofConf = 1;
00945 
00946   LOG_INFO << "[StBTofGeometry] Done. NValidTrays = " << mNValidTrays << "  NModulesInTray = " << mModulesInTray << endm;
00947 
00948   if(mDebug) Print();
00949 
00950   return;
00951 
00952 }
00953 
00954 //_____________________________________________________________________________
00955 Bool_t StBTofGeometry::ContainOthers(TVolume *element)
00956 {
00957   TVolumeView *elementView = new TVolumeView(*element,1);
00958   TList *list = elementView->GetListOfShapes();
00959   if ( list ) {
00960     LOG_INFO << " yes list in " << element->GetName() << endm;
00961     return kTRUE;
00962   } else {
00963     LOG_INFO << " no list in " << element->GetName() << endm;
00964     return kFALSE;
00965   }
00966 }
00967 
00968 //_____________________________________________________________________________
00969 Bool_t StBTofGeometry::LackThis(const char* fromWhere)
00970 {
00971    if (gBTofGeometry == 0) {
00972       LOG_INFO << " !! Warning from " << fromWhere
00973            << "\n no StBTofGeometry existing, create one instance first"
00974            << endm;
00975       return kTRUE;
00976    } else return kFALSE;
00977 }
00978 
00979 //_____________________________________________________________________________
00980 Int_t StBTofGeometry::CalcCellId(const Int_t volumeId, const Float_t* local)
00981 const
00982 {
00983    Double_t dlocal[3];
00984    dlocal[0] = local[0];
00985    dlocal[1] = local[1];
00986    dlocal[2] = local[2];
00987    return CalcCellId(volumeId,dlocal);
00988 }
00989 
00990 //_____________________________________________________________________________
00991 Int_t StBTofGeometry::CalcCellId(const Int_t volumeId, const Double_t* local)
00992 const
00993 {
00994    //
00995    // Calculate cellID based on volumeId and local position in a sensor
00996    //
00997 
00998    Int_t icell, imodule, itray;
00999    DecodeVolumeId(volumeId, imodule, itray);
01000    StBTofGeomSensor *sensor = GetGeomSensor(imodule, itray);
01001    if (!sensor) return -1;
01002    icell = sensor->FindCellIndex(local);
01003    Int_t ret = CalcCellId(icell, imodule, itray);
01004 
01005    return ret;
01006 }
01007 
01008 //_____________________________________________________________________________
01009 Bool_t StBTofGeometry::IsCellValid(const Int_t icell)
01010 const
01011 {
01012    //
01013    //Check the validity of cell# = [1,mCellsInModule]
01014    //
01015    return (icell>=1 && icell<=mCellsInModule);
01016 }
01017 
01018 //_____________________________________________________________________________
01019 Bool_t StBTofGeometry::IsSensorValid(const Int_t imodule)
01020 const
01021 {
01022    //
01023    //Check the validity of module# = [1,mModulesInTray]
01024    //
01025    return (imodule>=1 && imodule<=mModulesInTray);
01026 }
01027 
01028 //_____________________________________________________________________________
01029 Bool_t StBTofGeometry::IsTrayValid(const Int_t itray)
01030 const
01031 {
01032    //
01033    //Check the validity of tray#
01034    // return kTRUE if found in the tray list
01035    //
01036 
01037    Bool_t ret = kFALSE;
01038    for(int i=0;i<mNValidTrays;i++) {
01039      if(!mBTofTray[i]) continue;
01040      if(mBTofTray[i]->Index() == itray) {
01041        ret = kTRUE;
01042        break;
01043      }
01044    }
01045 
01046    return ret;
01047 }
01048 
01049 //_____________________________________________________________________________
01050 Int_t StBTofGeometry::CalcSensorId(const Int_t imodule, const Int_t itray)
01051 const
01052 {
01053    //
01054    //Encode imodule and itray into SensorId
01055    //
01056 
01057    Int_t sensorId = -1;
01058    if (!IsSensorValid(imodule)) return sensorId;
01059 
01060    Int_t idx = GetAtOfTray(itray);
01061 
01062    if (idx<0) return sensorId;
01063 
01064    sensorId = imodule-1 + mModulesInTray*idx;
01065 
01066    return sensorId;
01067 }
01068 
01069 //_____________________________________________________________________________
01070 Int_t StBTofGeometry::PrevCellId(const Int_t cellId)
01071 const
01072 {
01073    //
01074    //Look up the cell prior to cellId in the same module
01075    // and return the cellId of cell found
01076    //
01077 
01078    Int_t found = -1;
01079    Int_t icell, imodule, itray;
01080    DecodeCellId(cellId,icell,imodule,itray);
01081    StBTofGeomSensor* sensor = GetGeomSensor(imodule,itray);
01082    Int_t icell_p = sensor->PrevCellIndex(icell);
01083 
01084    found = CalcCellId(icell_p,imodule,itray);
01085 
01086    return found;
01087 }
01088 
01089 //_____________________________________________________________________________
01090 Int_t StBTofGeometry::NextCellId(const Int_t cellId)
01091 const
01092 {
01093    //
01094    //Look up the cell prior to cellId in the same module
01095    // and return the cellId of cell found
01096    //
01097 
01098    Int_t found = -1;
01099    Int_t icell, imodule, itray;
01100    DecodeCellId(cellId,icell,imodule,itray);
01101    StBTofGeomSensor* sensor = GetGeomSensor(imodule,itray);
01102    Int_t icell_p = sensor->NextCellIndex(icell);
01103 
01104    found = CalcCellId(icell_p,imodule,itray);
01105 
01106    return found;
01107 }
01108 
01109 //_____________________________________________________________________________
01110 Int_t StBTofGeometry::CalcCellId(const Int_t icell, const Int_t imodule,
01111                                  const Int_t itray)
01112 const
01113 {
01114    //
01115    //Encode icell, imodule and itray into CellId
01116    //
01117 
01118    Int_t cellId = -1;
01119    if (!IsCellValid(icell)) return cellId;
01120 
01121    Int_t sensorId = CalcSensorId(imodule,itray);
01122    if (sensorId<0) return cellId;                 //Invalid sensorId
01123    
01124    cellId = icell-1 + mCellsInModule*sensorId;
01125 
01126    return cellId;
01127 }
01128 
01129 //_____________________________________________________________________________
01130 Bool_t  StBTofGeometry::DecodeSensorId(const Int_t sensorId,
01131                                              Int_t &imodule, Int_t &itray)
01132 const
01133 {
01134    //
01135    //decode sensorId into tray#, module# in itray, imodule
01136    //
01137 
01138    imodule   = sensorId%mModulesInTray + 1;
01139    if (!IsSensorValid(imodule)) return kFALSE;
01140 
01141    Int_t idx = sensorId/mModulesInTray;
01142 
01143    itray = GetTrayIndexAt(idx);
01144 
01145    StBTofGeomTray *tray = GetGeomTrayAt(idx);
01146    if (!tray) return kFALSE;
01147    itray   = tray->Index();
01148 //   delete tray;
01149    return kTRUE;
01150 }
01151 
01152 //_____________________________________________________________________________
01153 Bool_t  StBTofGeometry::DecodeCellId(const Int_t cellId, Int_t &icell,
01154                                            Int_t &imodule, Int_t &itray)
01155 const
01156 {
01157    //
01158    //decode cellId into tray#, module#, cell# in itray, imodule, icell
01159    //
01160 
01161    Int_t sensorId = cellId/mCellsInModule;
01162    if (!DecodeSensorId(sensorId,imodule,itray)) return kFALSE;
01163 
01164    icell  = cellId%mCellsInModule + 1;
01165 
01166    return IsCellValid(icell);
01167 }
01168 
01169 //_____________________________________________________________________________
01170 Int_t StBTofGeometry::GetCellIndex(const Int_t cellId)
01171 const
01172 {
01173    //
01174    //decode cellId and return cell#
01175    //
01176 
01177    Int_t icell  = cellId%mCellsInModule + 1;
01178 
01179    return icell;
01180 }
01181 
01182 //_____________________________________________________________________________
01183 void  StBTofGeometry::DecodeVolumeId(const Int_t volumeId,
01184                                            Int_t &imodule, Int_t &itray)
01185 const
01186 {
01187    //
01188    //decode the volumeId into tray# and module# in itray, imodule
01189    // see the definition of TOFr's volumeId in "pams/sim/g2t/g2t_volume_id.g"
01190    //
01191 
01192    Int_t ires    = volumeId;
01193 
01194    Int_t rileft  = int(ires/10/100/100);
01195    ires          = ires-rileft*100*100*10;
01196 
01197    itray         = int(ires/10/100);
01198    ires          = ires-itray*100*10;
01199 
01200    imodule       = int(ires/10);
01201    // Int_t icell   = ires%10;  //always 0 since not defined in geant geometry
01202 
01203    itray = itray + (rileft-1)*mSectorsInBTOH;
01204 
01205    return;
01206 }
01207 
01208 //_____________________________________________________________________________
01209 StBTofGeomSensor* StBTofGeometry::GetGeomCell(const Int_t cellId)
01210 const
01211 {
01212    //
01213    //Return StBTofGeomSensor* where the cell of cellId is in
01214    //
01215 
01216    Int_t icell, imodule, itray;
01217    DecodeCellId(cellId, icell, imodule, itray);
01218    StBTofGeomSensor* sensor = GetGeomSensor(imodule, itray);
01219 
01220    return sensor;
01221 }
01222 
01223 //_____________________________________________________________________________
01224 StBTofGeomSensor* StBTofGeometry::GetGeomSensor(const Int_t imodule,
01225                                                 const Int_t itray)
01226 const
01227 {
01228    //
01229    //itray is dummy if itray==0 and it is the current single tray
01230    //
01231 
01232    for(int i=0;i<mNValidTrays;i++) {
01233      if(!mBTofTray[i]) continue;
01234      if(mBTofTray[i]->Index()==itray) {
01235        for(int j=0;j<mModulesInTray;j++) {
01236          if(!mBTofSensor[i][j]) continue;
01237          if(mBTofSensor[i][j]->Index()==imodule) return mBTofSensor[i][j];
01238        }
01239      }
01240    }
01241 
01242    return 0;
01243 }
01244 
01245 //_____________________________________________________________________________
01246 StBTofGeomTray*   StBTofGeometry::GetGeomTray(const Int_t itray)
01247 const
01248 {
01249    //
01250    //itray is dummy if itray==0 and it is the current single tray
01251    //
01252 
01253    for(int i=0;i<mNValidTrays;i++) {
01254      if(!mBTofTray[i]) continue;
01255      if(mBTofTray[i]->Index()==itray) return mBTofTray[i];
01256    }
01257    return 0;
01258 }
01259 
01260 //_____________________________________________________________________________
01261 Int_t StBTofGeometry::GetTrayIndexAt(const Int_t idx)
01262 const
01263 {
01264    //
01265    //Get the tray index at index of the list
01266    //
01267 
01268    TDataSetIter nextSector(mTopNode);
01269    TVolumeView *sectorVolume = 0;
01270    Int_t ibtoh = 0, i = 0;
01271    Int_t itray = -1;
01272    while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
01273      i++;
01274      if ( i>mSectorsInBTOH ) ibtoh = 1;
01275      int trayIndex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
01276      if (i==idx) {
01277        itray = trayIndex;
01278        break;
01279      }
01280    }
01281 
01282    return itray;
01283 }
01284 
01285 //_____________________________________________________________________________
01286 StBTofGeomTray* StBTofGeometry::GetGeomTrayAt(const Int_t idx)
01287 const
01288 {
01289    //
01290    //Get the StBTofGeomTray at index of the list
01291    //
01292 
01293    Int_t itray = GetTrayIndexAt(idx);
01294    if(itray<=0||itray>mNTrays) return 0;
01295 
01296    for(int i=0;i<mNValidTrays; i++) {
01297      if(mBTofTray[i] && mBTofTray[i]->Index()==itray) return mBTofTray[i];
01298    }
01299 
01300    return 0;
01301 }
01302 
01303 //_____________________________________________________________________________
01304 Int_t  StBTofGeometry::GetAtOfTray(const Int_t itray)
01305 const
01306 {
01307    //
01308    //Find out the list-index of StBTofGeomTray with TrayIndex=itray
01309    // itray is dummy if itray==0 and it is the current single tray
01310    //
01311 
01312    Int_t at = -1;
01313 
01314    TDataSetIter nextSector(mTopNode);
01315    TVolumeView *sectorVolume = 0;
01316    Int_t ibtoh = 0, i = 0;
01317    while ( (sectorVolume = (TVolumeView *)nextSector())  ) {
01318      i++;
01319      if ( i>mSectorsInBTOH ) ibtoh = 1;
01320      int trayIndex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
01321      if (trayIndex == itray) {
01322        at = i;
01323        break;
01324      }
01325    }
01326 
01327    return at;
01328 }
01329 
01330 //_____________________________________________________________________________
01331 void StBTofGeometry::Print(Option_t *opt)  const
01332 {
01333    LOG_INFO << "Trays=" << mNValidTrays <<"\t ModulesInTray=" << mModulesInTray
01334         << "\t CellsInModule=" << mCellsInModule << endm;
01335 }
01336 
01337 //_____________________________________________________________________________
01338 Int_t StBTofGeometry::CellIdPointIn(const StThreeVectorD& point)
01339 const
01340 {
01341    //
01342    //Return the ID of cell which containing the global point
01343    //
01344 
01345    Int_t cellId = -1;
01346    Double_t xl[3], xg[3];
01347    xg[0] = point.x();
01348    xg[1] = point.y();
01349    xg[2] = point.z();
01350 
01351    //Loop over trays
01352    //
01353 
01354    Int_t itray = -1, imodule = -1, icell = -1;
01355    for(int i=0;i<mNValidTrays;i++) {
01356      if(!mBTofTray[i]) continue;
01357      if ( mBTofTray[i]->IsGlobalPointIn(point) ) {
01358        itray = mBTofTray[i]->Index();
01359        if ( !(mBTofTray[i]->GetfView()->GetListSize()) ) {
01360          LOG_INFO << " No sensors in tray " << itray << endm;
01361          return cellId;
01362        }
01363        
01364        for( int j=0;j<mModulesInTray;j++) {
01365          if(!mBTofSensor[i][j]) continue;
01366          if ( mBTofSensor[i][j]->IsGlobalPointIn(point) ) {
01367            imodule = mBTofSensor[i][j]->Index();
01368            mBTofSensor[i][j]->Master2Local(xg,xl);
01369            icell = mBTofSensor[i][j]->FindCellIndex(xl);
01370          }
01371        } // end for (j)
01372      } // end if
01373    } // end for (i)
01374 
01375    if ( itray <= 0 || imodule <= 0 ) return cellId;
01376    cellId = CalcCellId(icell, imodule, itray);
01377    return cellId;
01378 }
01379 
01380 //_____________________________________________________________________________
01381 Bool_t StBTofGeometry::HelixCross(const StHelixD &helix)
01382 const
01383 {
01384    //
01385    // return "kTRUE" if any cell is crossed by this helix
01386    //
01387 
01388    IntVec idVec;
01389    DoubleVec pathVec;
01390    PointVec crossVec;
01391 
01392    Bool_t crossed = HelixCrossCellIds(helix,idVec,pathVec,crossVec);
01393 
01394    return crossed;
01395 }
01396 
01397 //_____________________________________________________________________________
01398 Bool_t StBTofGeometry::HelixCrossCellIds(const StHelixD &helix,
01399                        IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec)
01400 const
01401 {
01402    //
01403    // return "kTRUE" if any cell is crossed by this helix
01404    //  and also fill the cellIds which are crossed by the helix
01405    //  and the path length of the helix before crossing the cell
01406    //
01407 
01408    IntVec projTrayVec;
01409    if( !projTrayVector(helix, projTrayVec) ) return kFALSE;
01410 
01411    Double_t pathLen;
01412    Int_t cellId;
01413    StThreeVectorD cross;
01414    idVec.clear();
01415    pathVec.clear();
01416    crossVec.clear();
01417 
01418    for(int i=0;i<mNValidTrays;i++) {
01419      if(!mBTofTray[i]) continue;
01420      int trayId = mBTofTray[i]->Index();
01421 
01423 
01424      bool itrayFind = kFALSE;
01425      for(size_t it=0;it<projTrayVec.size();it++) {
01426        int validtrayId = projTrayVec[it];
01427        if(validtrayId==trayId) {
01428          itrayFind = kTRUE;
01429          break;
01430        }
01431      }
01432      if(!itrayFind) continue;
01433 
01434      for(int j=0;j<mModulesInTray;j++) {
01435        if(!mBTofSensor[i][j]) continue;
01436        int moduleId = mBTofSensor[i][j]->Index();
01437        if ( mBTofSensor[i][j]->HelixCross(helix,pathLen,cross) ) {         
01438          Double_t global[3], local[3];
01439          global[0] = cross.x();
01440          global[1] = cross.y();
01441          global[2] = cross.z();
01442          mBTofSensor[i][j]->Master2Local(global,local);
01443          Int_t icell = mBTofSensor[i][j]->FindCellIndex(local);
01444          cellId = CalcCellId(icell, moduleId, trayId);
01445          if (cellId>=0) {    // a bug before // reject hit in the edge of module;
01446            pathVec.push_back(pathLen);
01447            idVec.push_back(cellId);
01448            crossVec.push_back(cross);
01449          }
01450        }
01451      } // end for (j)
01452    } // end for (i)
01453    
01454    if (idVec.size()>0) {
01455      //     mleak1.PrintMem(" normal return true");
01456      return kTRUE;
01457    }
01458    else {
01459      //     mleak1.PrintMem(" normal return false");
01460      return kFALSE;
01461    }
01462 }
01463 
01464 //_____________________________________________________________________________
01465 Bool_t StBTofGeometry::HelixCrossCellIds(const StHelixD &helix,
01466                        IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec, DoubleVec &thetaVec)
01467 const
01468 {
01469    //
01470    // return "kTRUE" if any cell is crossed by this helix
01471    //  and also fill the cellIds which are crossed by the helix
01472    //  and the path length of the helix before crossing the cell
01473    //
01474 
01475    IntVec projTrayVec;
01476    if( !projTrayVector(helix, projTrayVec) ) return kFALSE;
01477 
01478    Double_t pathLen,theta;
01479    Int_t cellId;
01480    StThreeVectorD cross;
01481    idVec.clear();
01482    pathVec.clear();
01483    crossVec.clear();
01484 
01485    for(int i=0;i<mNValidTrays;i++) {
01486      if(!mBTofTray[i]) continue;
01487      int trayId = mBTofTray[i]->Index();
01488 
01490 
01491      bool itrayFind = kFALSE;
01492      for(size_t it=0;it<projTrayVec.size();it++) {
01493        int validtrayId = projTrayVec[it];
01494        if(validtrayId==trayId) {
01495          itrayFind = kTRUE;
01496          break;
01497        }
01498      }
01499      if(!itrayFind) continue;
01500 
01501      for(int j=0;j<mModulesInTray;j++) {
01502        if(!mBTofSensor[i][j]) continue;
01503        int moduleId = mBTofSensor[i][j]->Index();
01504        if ( mBTofSensor[i][j]->HelixCross(helix,pathLen,cross,theta) ) {
01505          Double_t global[3], local[3];
01506          global[0] = cross.x();
01507          global[1] = cross.y();
01508          global[2] = cross.z();
01509          mBTofSensor[i][j]->Master2Local(global,local);
01510          Int_t icell = mBTofSensor[i][j]->FindCellIndex(local);
01511          cellId = CalcCellId(icell, moduleId, trayId);
01512          if (cellId>=0) {    // a bug before // reject hit in the edge of module;
01513            pathVec.push_back(pathLen);
01514            idVec.push_back(cellId);
01515            crossVec.push_back(cross);
01516            thetaVec.push_back(theta);
01517 
01518          }
01519        }
01520      } // end for (j)
01521    } // end for (i)
01522 
01523    if (idVec.size()>0) {
01524      //     mleak1.PrintMem(" normal return true");
01525      return kTRUE;
01526    }
01527    else {
01528      //     mleak1.PrintMem(" normal return false");
01529      return kFALSE;
01530    }
01531 }
01532 
01533 //_____________________________________________________________________________
01534 Bool_t StBTofGeometry::HelixCross(const StHelixD &helix, IntVec validModuleVec, IntVec projTrayVec)
01535 const
01536 {
01537    //
01538    // return "kTRUE" if any cell is crossed by this helix
01539    //
01540 
01541    IntVec idVec;
01542    DoubleVec pathVec;
01543    PointVec crossVec;
01544 
01545    Bool_t crossed = HelixCrossCellIds(helix,validModuleVec, projTrayVec, idVec,pathVec,crossVec);
01546 
01547    return crossed;
01548 }
01549 
01550 
01551 //_____________________________________________________________________________
01552 Bool_t StBTofGeometry::HelixCrossCellIds(const StHelixD &helix, IntVec validModuleVec, IntVec projTrayVec, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec)
01553 const
01554 {
01556   // optimized :
01557   // input : helix, validModuleVec from DAQ, 
01558   //         projTrayVec possible hit tray from this helix
01559   //                                   Xin Dong
01561    //
01562    // return "kTRUE" if any cell is crossed by this helix
01563    //  and also fill the cellIds which are crossed by the helix
01564    //  and the path length of the helix before crossing the cell
01565    //
01566 
01567    if(validModuleVec.size()==0) return kFALSE;
01568    if(projTrayVec.size()==0) return kFALSE;
01569 
01570    Double_t pathLen;
01571    Int_t cellId;
01572    StThreeVectorD cross;
01573    idVec.clear();
01574    pathVec.clear();
01575    crossVec.clear();
01576 
01577    for(int i=0;i<mNValidTrays;i++) {
01578      if(!mBTofTray[i]) continue;
01579      int trayId = mBTofTray[i]->Index();
01580      bool itrayFind = kFALSE;
01581 
01582      for(size_t it=0;it<projTrayVec.size();it++) {
01583        int validtrayId = projTrayVec[it];
01584        if(validtrayId==trayId) {
01585          itrayFind = kTRUE;
01586          break;
01587        }
01588      }
01589      if(!itrayFind) continue;
01590 
01591      //     LOG_INFO << " Helix cross sensitive tray " << trayId << endm;
01592 
01593      for(int j=0;j<mModulesInTray;j++) {
01594        if(!mBTofSensor[i][j]) continue;
01595        int moduleId = mBTofSensor[i][j]->Index();
01596        for(size_t iv=0;iv<validModuleVec.size();iv++) {
01597          int validtrayId = validModuleVec[iv]/100;
01598          int validmoduleId = validModuleVec[iv]%100;
01599          if(validtrayId==trayId&&validmoduleId==moduleId) {
01600            if ( mBTofSensor[i][j]->HelixCross(helix,pathLen,cross) ) {     
01601              Double_t global[3], local[3];
01602              global[0] = cross.x();
01603              global[1] = cross.y();
01604              global[2] = cross.z();
01605              mBTofSensor[i][j]->Master2Local(global,local);
01606              Int_t icell = mBTofSensor[i][j]->FindCellIndex(local);
01607              cellId = CalcCellId(icell, moduleId, trayId);
01608              if (cellId>=0) {   // a bug before // reject hit in the edge of module;
01609                pathVec.push_back(pathLen);
01610                idVec.push_back(cellId);
01611                crossVec.push_back(cross);
01612              }
01613            }
01614          } // endif (tray && module)
01615        } // end for (iv)
01616      } // end for (j)
01617    } // end for (i)
01618 
01619    if (idVec.size()>0) {
01620      //     mleak1.PrintMem(" normal return true");
01621      return kTRUE;
01622    }
01623    else {
01624      //     mleak1.PrintMem(" normal return false");
01625      return kFALSE;
01626    }
01627 }
01628 
01629 //---------------------------------------------------------------------------
01630 // estimate the possible projection on the TOF tray
01631 Bool_t StBTofGeometry::projTrayVector(const StHelixD &helix, IntVec &trayVec) const {
01632 
01633   trayVec.clear();
01634   double R_tof[2]= {210., 216.};  // inner and outer surfaces
01635 //  double res = 5.0;
01636 // second choice - choose 2 neighbouring trays
01637 
01638   for(int i=0;i<2;i++) {
01639 
01640     double s = helix.pathLength(R_tof[i]).first;
01641     if(s<0.) s = helix.pathLength(R_tof[i]).second;
01642     StThreeVectorD point = helix.at(s);
01643     double phi = point.phi()*180/3.14159;
01644     double z = point.z();
01645 
01646     int itray[3] = {0,0,0};
01647     if(z<0) {
01648       // east ring, start from 108 deg (id=61) , clock-wise from east facing west
01649       itray[0] = (255+(int)phi)%360/6+61;
01650       itray[1] = itray[0] - 1;
01651       if(itray[1]<=60) itray[1] += 60;
01652       itray[2] = itray[0] + 1;
01653       if(itray[2]>120) itray[2] -= 60;
01654 //      itray[1] = (255+(int)(phi+res))%360/6+61;
01655 //      itray[2] = (255+(int)(phi-res))%360/6+61;
01656     } else {
01657       // west ring, start from 72 deg (id=1) , clock-wise from west facing east
01658       itray[0] = (435-(int)phi)%360/6+1;
01659       itray[1] = itray[0] - 1;
01660       if(itray[1]<=0) itray[1] += 60;
01661       itray[2] = itray[0] + 1;
01662       if(itray[2]>60) itray[2] -= 60;
01663 //      itray[1] = (435-(int)(phi+res))%360/6+1;
01664 //      itray[2] = (435-(int)(phi-res))%360/6+1;
01665     }
01666 
01667     for(int k=0;k<3;k++) {
01668       if(itray[k]<=0 || itray[k]>120) continue;
01669 
01670       bool found = kFALSE;
01671       for(size_t j=0;j<trayVec.size();j++) {
01672         if(trayVec[j]==itray[k]) {
01673           found = kTRUE;
01674           break;
01675         }
01676       }
01677 
01678       if(found) continue;
01679 
01680       trayVec.push_back(itray[k]);
01681     } // end loop k  
01682   } // end loop i
01683 
01684 /*
01685    cout << " proj tray id = ";
01686    for(size_t it=0;it<trayVec.size();it++) {
01687      cout << trayVec[it] << " ";
01688    }
01689    cout << endl;
01690 */
01691   
01692   if(trayVec.size()>0) return kTRUE;
01693   else return kFALSE;
01694 }