loon: GeoShieldGroup.cxx Source File

00001 
00002 //
00003 // GeoShieldGroup
00004 //
00005 // GeoShieldGroup is a helper class for managing a single shield group
00006 //
00007 // Author:  S. Kasahara 10/05
00009 
00010 #include <iostream>
00011 using namespace std;
00012 
00013 #include "TGeoVolume.h"
00014 #include "TGeoMatrix.h"
00015 #include "TGeoXtru.h"
00016 #include "TGeoBBox.h"
00017 #include "GeoGeometry/GeoGeometry.h"
00018 #include "GeoGeometry/GeoShieldGroup.h"
00019 #include "GeoGeometry/GeoMediumMap.h"
00020 #include "GeoGeometry/GeoMedium.h"
00021 #include "Plex/PlexPlaneId.h"
00022 #include "UgliGeometry/UgliDbiSteelPln.h"
00023 #include "MessageService/MsgService.h"
00024 
00025 ClassImp(GeoShieldGroup)
00026 
00027 CVSID("$Id: GeoShieldGroup.cxx,v 1.8 2008/01/09 06:38:11 schubert Exp $");
00028 
00029 //_____________________________________________________________________________
00030 GeoShieldGroup::GeoShieldGroup(GeoGeometry* geo, 
00031                                GeoShield::EGroupType group): 
00032                                fGeoGeometry(geo), fGroup(group),
00033                                fPoints(),fVolumes() {
00034   // Normal constructor
00035 
00036   MSG("Geo",Msg::kDebug) << "GeoShieldGroup normal ctor @ " << this << endl;
00037   for ( int ic = 0; ic < 3; ic++ ) {
00038     fPoints[0][ic] = +999999.;
00039     fPoints[1][ic] = -999999.;
00040   }
00041   
00042 }
00043 
00044 //_____________________________________________________________________________
00045 GeoShieldGroup::~GeoShieldGroup() {
00046   // Destructor, delete all owned objects
00047   MSG("Geo",Msg::kDebug) << "GeoShield dtor @ " << this << endl;
00048   
00049   // geogeometry creator is not owned
00050   // volumes in list are not owned
00051 
00052 }
00053 
00054 //_____________________________________________________________________________
00055 void GeoShieldGroup::Print(Option_t* /* option */) const {
00056    // print to cout
00057 
00058   if ( fVolumes.size() == 0 ) {
00059     cout << "Group " << GeoShield::AsString(fGroup) << " has no volumes." 
00060          << endl;
00061   }
00062   else {
00063     cout << "Group " << GeoShield::AsString(fGroup) << " has " 
00064          << fVolumes.size() << " volume(s):" << endl;
00065     std::list<const TGeoVolume*>::const_iterator volitr;
00066     for ( volitr = fVolumes.begin(); volitr != fVolumes.end(); volitr++ ) {
00067        (*volitr) -> Print();
00068     }
00069   }
00070   
00071   cout << "Group box boundary points in global coordinates: " << endl;
00072   for ( int ipt = 0; ipt < 2; ipt++ ) {
00073     cout << " " << ipt << ": (" 
00074          << fPoints[ipt][0] << "," << fPoints[ipt][1] 
00075          << "," << fPoints[ipt][2] << ")" << endl;
00076   }
00077   
00078 }
00079 
00080 //_____________________________________________________________________________
00081 void GeoShieldGroup::AddVolume(const TGeoVolume* pairVol, 
00082                                const UgliDbiSteelPln* stRow) {
00083   // The stRow is used to extract positional information in order to
00084   // determine size of bounding box
00085 
00086   fVolumes.push_back(pairVol);
00087 
00088   Double_t scale = Geo::GetScale(fGeoGeometry->GetAppType());
00089   
00090   // Extract size of volume 
00091   // Convert to local coordinates of pair to add dimensions of volume,
00092   // then convert these coordinates to global to determine size of
00093   // group, stored in global coordinates
00094   Double_t gxyz[3] = {(stRow->GetX0())*scale,(stRow->GetY0())*scale,
00095                       (stRow->GetZBack()-0.5*(stRow->GetTotalZ()))*scale};
00096   TGeoCombiTrans* pairMatrix = new TGeoCombiTrans(pairVol->GetName(),
00097                                    gxyz[0],gxyz[1],gxyz[2],new TGeoRotation(
00098                                    pairVol->GetName(),
00099                                    stRow->GetThetaDeg(0),stRow->GetPhiDeg(0),
00100                                    stRow->GetThetaDeg(1),stRow->GetPhiDeg(1),
00101                                    stRow->GetThetaDeg(2),stRow->GetPhiDeg(2)));
00102   pairMatrix -> RegisterYourself();
00103 
00104   Double_t lxyz[3] = {0};  
00105   pairMatrix -> MasterToLocal(gxyz,lxyz);
00106 
00107   // Use the volume bbox to determine dimensions of volume
00108   TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(pairVol -> GetShape());
00109   Double_t halfw[3] = {bbox->GetDX(),bbox->GetDY(),bbox->GetDZ()};
00110   Double_t pad = 0.; // 0.05*scale; // not necessary
00111 
00112   // In local group coordinates, store the point, storage depends on
00113   // group type
00114   switch ( fGroup ) {
00115 
00116   case GeoShield::kFarInnerE:
00117   case GeoShield::kFarInnerW:
00118     {
00119       Double_t lxyzmin[3] = {0};
00120       Double_t lxyzmax[3] = {0};
00121       Double_t gxyzmin[3] = {0};
00122       Double_t gxyzmax[3] = {0};
00123       for ( int ic = 0; ic < 3; ic++ ) {
00124         lxyzmin[ic] = lxyz[ic] - halfw[ic] - pad;
00125         lxyzmax[ic] = lxyz[ic] + halfw[ic] + pad;
00126         pairMatrix -> LocalToMaster(lxyzmin,gxyzmin);
00127         pairMatrix -> LocalToMaster(lxyzmax,gxyzmax);
00128         fPoints[0][ic] = TMath::Min(fPoints[0][ic],gxyzmin[ic]);
00129         fPoints[0][ic] = TMath::Min(fPoints[0][ic],gxyzmax[ic]);
00130         fPoints[1][ic] = TMath::Max(fPoints[1][ic],gxyzmin[ic]);
00131         fPoints[1][ic] = TMath::Max(fPoints[1][ic],gxyzmax[ic]);
00132       }
00133     }
00134   case GeoShield::kFarOuterE: 
00135   case GeoShield::kFarOuterW: 
00136   case GeoShield::kFarTop:
00137     {
00138       Double_t lxyzcnr[3] = {0};
00139       Double_t gxyzcnr[3] = {0};
00140       Double_t sign[2] = {-1,1};
00141       for ( int ip = 0; ip < 8; ip++ ) {
00142         for ( int ix = 0; ix < 2; ix++ ) {
00143           lxyzcnr[0] = lxyz[0] + sign[ix]*(halfw[0] + pad);
00144           for ( int iy = 0; iy < 2; iy++ ) {
00145             lxyzcnr[1] = lxyz[1] + sign[iy]*(halfw[1] + pad);
00146             for ( int iz = 0; iz < 2; iz++ ) {
00147               lxyzcnr[2] = lxyz[2] + sign[iz]*(halfw[2] + pad);   
00148               pairMatrix -> LocalToMaster(lxyzcnr,gxyzcnr);
00149               for ( int ic = 0; ic < 3; ic++ ) {
00150                 fPoints[0][ic] = TMath::Min(fPoints[0][ic],gxyzcnr[ic]);
00151                 fPoints[1][ic] = TMath::Max(fPoints[1][ic],gxyzcnr[ic]);
00152               }
00153             }
00154           }
00155         }
00156       }
00157     }
00158     break;
00159 
00160   default:
00161     MSG("Geo",Msg::kFatal) << "ShieldGroup " 
00162                            <<  GeoShield::AsString(fGroup)
00163                            << " unknown!. Abort." << endl;
00164     abort();
00165   } // end of shield group switch
00166 
00167   return;
00168   
00169 }
00170 
00171 //_____________________________________________________________________________
00172 void GeoShieldGroup::BuildNode(TGeoVolume* hallVol) {
00173   // Build shield group TGeoNode and add to the hall volume.
00174   // All shield volumes are added to the shield group volume as
00175   // nodes as well.
00176 
00177   Double_t scale = Geo::GetScale(fGeoGeometry->GetAppType());
00178 
00179   // All shield group bounding boxes are made up of air, i.e.
00180   // blend into hall mother volume
00181   const GeoMediumMap& medMap = fGeoGeometry->GetMediumMap();
00182  
00183   // Shield group volume positions are created relative to MARS, so
00184   // need hall global position to recalculate.
00185   TGeoNode* linrNode = gGeoManager->GetVolume("MARS")->GetNode("LINR_1");
00186   const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
00187   
00188   // So far, all shield group volumes are built as TGeoBBox
00189   Double_t gxyz0[3] = {0}; // global box center
00190   for ( int ip = 0; ip < 2; ip++ ) {
00191     for ( int ic = 0; ic < 3; ic++ ) {
00192       gxyz0[ic] += fPoints[ip][ic];
00193     }
00194   }
00195   for ( int ic = 0; ic < 3; ic++ ) gxyz0[ic] /= 2.;
00196 
00197   Double_t halfxyz[3] = {0}; // box halfwidths
00198   for ( int ip = 0; ip < 2; ip++ ) {
00199     for ( int ic = 0; ic < 3; ic++ ) {
00200       halfxyz[ic] 
00201       = TMath::Max(halfxyz[ic],TMath::Abs(fPoints[ip][ic]-gxyz0[ic]));
00202     }
00203   }
00204 
00205   // Make shield group box volume
00206   std::string grpName = GeoShield::AsString(fGroup);
00207   TGeoVolume* grpVol = 0;
00208   
00209   switch ( fGroup ) {
00210 
00211   case GeoShield::kFarInnerE:
00212   case GeoShield::kFarInnerW:
00213   case GeoShield::kFarOuterE:
00214   case GeoShield::kFarOuterW:
00215   {     
00216     // inner and outer wings are simple boxes
00217     grpVol = gGeoManager->MakeBox(grpName.c_str(),medMap.GetMedium(Geo::kHall),
00218                           halfxyz[0],halfxyz[1],halfxyz[2]);
00219   }
00220   break;
00221   
00222   case GeoShield::kFarTop:
00223   {
00224     // top section is a tgeoxtru outline of a box with trapezoid cut-out to
00225     // fit over top of detector
00226     TGeoXtru* grpShp = new TGeoXtru(2); // odd not to have a ctor w/name arg
00227     grpShp -> SetName(grpName.c_str());
00228 
00229     // This little bit of math is affordable because it's only done once
00230     Double_t hiy = 4.25*scale;
00231     Double_t loy = gxyz0[1] - halfxyz[1];
00232     Double_t tand225 = 0.414213562;
00233     Double_t halfxhiy = hiy*tand225;
00234     Double_t halfxloy = (hiy-loy)+halfxhiy;
00235     Double_t halfy = (hiy-loy)/2.;
00236 
00237     const int nv = 8;
00238     Double_t x[nv] = {0};
00239     Double_t y[nv] = {0};
00240     
00241     // etch'a sketch
00242     x[0] = -halfxyz[0]; y[0] = -halfxyz[1];
00243     x[1] = -halfxyz[0]; y[1] = +halfxyz[1];
00244     x[2] = +halfxyz[0]; y[2] = +halfxyz[1];
00245     x[3] = +halfxyz[0]; y[3] = -halfxyz[1];
00246     x[4] = +halfxloy;   y[4] = -halfxyz[1];
00247     x[5] = +halfxhiy;   y[5] = -halfxyz[1]+2.*halfy;
00248     x[6] = -halfxhiy;   y[6] = -halfxyz[1]+2.*halfy;
00249     x[7] = -halfxloy;   y[7] = -halfxyz[1];
00250         
00251     grpShp -> DefinePolygon(8,x,y);
00252     // Define section must be called after DefinePolygon
00253     grpShp -> DefineSection(0,-halfxyz[2]);
00254     grpShp -> DefineSection(1,+halfxyz[2]);
00255     
00256     grpVol = new TGeoVolume(grpName.c_str(),grpShp,
00257                             medMap.GetMedium(Geo::kHall));  
00258   }
00259   break;
00260   
00261   default:
00262     MSG("Geo",Msg::kFatal) << "ShieldGroup " 
00263                            <<  GeoShield::AsString(fGroup)
00264                            << " unknown!. Abort." << endl;
00265     abort();
00266   } // end of shield group switch
00267 
00268   //grpVol -> SetVisibility(kTRUE);
00269   //grpVol -> SetLineColor(kBlack);
00270   grpVol -> SetVisibility(kFALSE);
00271   grpVol -> VisibleDaughters(kTRUE);
00272   
00273   // Place in hallVol as node
00274   // TGeoTranslation is adopted by gGeoManager
00275   // overlap because overlaps pair plane boxes. Fix me.
00276   hallVol -> AddNodeOverlap(grpVol,1,new TGeoTranslation(grpName.c_str(),
00277               (gxyz0[0]-h0[0]),(gxyz0[1]-h0[1]),(gxyz0[2]-h0[2])));
00278 
00279   // Now add volumes in group list to shield group volume
00280   std::list<const TGeoVolume*>::const_iterator volItr;
00281 
00282   for ( volItr = fVolumes.begin(); volItr != fVolumes.end(); volItr++ ) {
00283     const TGeoVolume* pairVol = *volItr;
00284     std::string pairName = pairVol -> GetName();
00285     TGeoCombiTrans* pairMatrix = dynamic_cast<TGeoCombiTrans*>
00286            (gGeoManager->GetListOfMatrices()->FindObject(pairName.c_str()));
00287     
00288     // Adjust pair matrix translation to set plane pair relative to volume
00289     const Double_t* x0pair = pairMatrix -> GetTranslation();
00290     pairMatrix -> SetTranslation(x0pair[0]-gxyz0[0],x0pair[1]-gxyz0[1],
00291                                  x0pair[2]-gxyz0[2]);
00292     
00293     grpVol -> AddNode(pairVol,1,pairMatrix);
00294   }  
00295   
00296   return;
00297   
00298 }
00299 
00300