Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]   [ 15.*.* ]  
001 #include "MuonGeoModel/BarrelToroidBuilderRDB.h"
002 #include "AthenaKernel/getMessageSvc.h"
003 
004 #include "RDBAccessSvc/IRDBRecord.h"
005 #include "RDBAccessSvc/IRDBRecordset.h"
006 #include "RDBAccessSvc/IRDBAccessSvc.h"
007                                                         
008 #include "MuonGeoModel/ArrayFunction.h"
009 
010 #include "GeoModelKernel/GeoBox.h"
011 #include "GeoModelKernel/GeoTube.h"
012 #include "GeoModelKernel/GeoTubs.h"
013 #include "GeoModelKernel/GeoPcon.h"
014 #include "GeoModelKernel/GeoCons.h"
015 #include "GeoModelKernel/GeoTrd.h"
016 #include "GeoModelKernel/GeoTrap.h"
017 #include "GeoModelKernel/GeoPara.h"
018 #include "GeoModelKernel/GeoPgon.h"
019 #include "GeoModelKernel/GeoMaterial.h"
020 #include "GeoModelKernel/GeoLogVol.h"
021 #include "GeoModelKernel/GeoPhysVol.h"
022 #include "GeoModelKernel/GeoFullPhysVol.h"
023 #include "GeoModelKernel/GeoTransform.h"
024 #include "GeoModelKernel/GeoAlignableTransform.h"
025 #include "GeoModelKernel/GeoNameTag.h"
026 #include "GeoModelKernel/GeoShapeShift.h"
027 #include "GeoModelKernel/GeoShapeUnion.h"
028 #include "GeoModelKernel/GeoShapeSubtraction.h"
029 #include "GeoModelKernel/GeoSerialTransformer.h" 
030 #include "GeoModelKernel/GeoIdentifierTag.h"
031 #include "GeoModelSvc/StoredMaterialManager.h"
032 
033 #include "StoreGate/StoreGateSvc.h"
034 
035 #include "CLHEP/GenericFunctions/Variable.hh"
036 
037 #include <stdexcept>
038 #include <vector>
039 #include <iomanip>
040 
041 #include <sstream>
042 typedef std::stringstream  my_sstream;
043 typedef std::ostringstream my_osstream;
044 
045 using namespace Genfun;
046 using namespace GeoXF;
047 
048 
049 namespace MuonGM {
050 BarrelToroidBuilderRDB::BarrelToroidBuilderRDB( StoreGateSvc  *pDetStore,
051                                                 IRDBAccessSvc *pRDBAccess, 
052                                                 std::string    geoTag,
053                                                 std::string    geoNode )    :
054     m_pRDBAccess(pRDBAccess), 
055     m_pDetStore (pDetStore)
056 {
057   MsgStream log( Athena::getMessageSvc(), "MuGM:BarrelToroBuildRDB" );
058 
059   log  <<  MSG::INFO  <<  "Fetching data with tag <"  <<  geoTag <<  ">"  <<  endreq;
060   m_Abrt = pRDBAccess->getRecordset("ABRT", geoTag, geoNode);
061   if (m_Abrt->size() == 0) {
062       log<<MSG::WARNING<<"Table ABRT not found in tag <"  <<  geoTag <<  ">"  <<" reading table ABRT-00" <<endreq;
063       m_Abrt = pRDBAccess->getRecordset("ABRT","ABRT-00");
064   }
065   m_Aect = pRDBAccess->getRecordset("AECT", geoTag, geoNode);
066   if (m_Aect->size() == 0) {
067       log<<MSG::WARNING<<"Table ABRT not found in tag <"  <<  geoTag <<  ">"  <<" reading table AECT-00" <<endreq;
068       m_Aect = pRDBAccess->getRecordset("AECT","AECT-00");
069   }
070   m_Feet = pRDBAccess->getRecordset("FEET", geoTag, geoNode);
071   if (m_Feet->size() == 0) {
072       log<<MSG::WARNING<<"Table FEET not found in tag <"  <<  geoTag <<  ">"  <<" reading table FEET-00" <<endreq;
073       m_Feet = pRDBAccess->getRecordset("FEET","FEET-00");
074   }
075   m_Rail = pRDBAccess->getRecordset("RAIL", geoTag, geoNode);
076   if (m_Rail->size() == 0) {
077       log<<MSG::WARNING<<"Table RAIL not found in tag <"  <<  geoTag <<  ">"  <<" reading table RAIL-00" <<endreq;
078       m_Rail = pRDBAccess->getRecordset("RAIL","RAIL-00");
079   }
080   
081   std::string Iron = "Iron";
082   std::string Aluminium = "Aluminium";
083 
084 
085 }
086 
087   void BarrelToroidBuilderRDB::buildBTVacuumVessel( GeoPhysVol* container ) 
088 {
089   //------------------------------------------------------------------------------------------
090   //  Builds BT Vacuum Vessel (tubes, ribs, and voussoir attachments)
091   //------------------------------------------------------------------------------------------
092 
093   const StoredMaterialManager*  theMaterialManager;
094   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
095   {
096     return;
097   } 
098 
099   const int maxDim = 8; 
100 
101 //---------
102    const double rMinBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMIN") * mm;                        
103    const double rMaxBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMAX") * mm;  
104    const double zMaxBTVessel      = (*m_Abrt)[0]->getFloat("CRYOZMAX") * mm;         
105    const double rCurvBTVessel     = (*m_Abrt)[0]->getFloat("CRYORCUR") * mm;                        
106    const double radiusBTVessel    = (*m_Abrt)[0]->getFloat("CRYORADI") * mm;       
107    const double thicknessBTVessel = (*m_Abrt)[0]->getFloat("CRYORADT") * mm;                        
108    const double zPosBTCryoring    = -1030.0 * mm; //(*m_Abrt)[0]->getFloat("CRYRNGZM") * mm;
109    const double radiusBTCryoring  = (*m_Abrt)[0]->getFloat("CRYRNGRA") * mm;
110                                                                               //
111    const std::string vesselMaterial = getMaterial("Iron"); 
112 //------------
113    int nBTRibs = 7;                                                  
114    double zPosBTRib[maxDim];                                                  //
115  
116  for (int i = 0; i<nBTRibs; i++)
117     zPosBTRib[i] = (*m_Abrt)[0]->getFloat("ZRIB",i) * mm;        
118                                                                        //
119    const double heightBTRib = (*m_Abrt)[0]->getFloat("CRYRIBYW") * mm;                              
120    const double zWidthBTRib = (*m_Abrt)[0]->getFloat("CRYRIBZL") * mm; 
121    const double CryRiWYp    = (*m_Abrt)[0]->getFloat("CRYRIWYP") * mm;
122    const double CryRiWYn    = (*m_Abrt)[0]->getFloat("CRYRIWYN") * mm;
123    const double CryRiWXp    = (*m_Abrt)[0]->getFloat("CRYRIWXP") * mm;
124    const double CryRiWXn    = (*m_Abrt)[0]->getFloat("CRYRIWXN") * mm;
125    const double CryRiWTh    = (*m_Abrt)[0]->getFloat("CRYRIWTH") * mm;
126                                                                                 
127 //----------------------------------------------------------------------------//
128 
129    const double DiamHoleOutLongTube   = (*m_Abrt)[0]->getFloat("CRYATTD0") * mm; //use also for building vosoir attachment
130    const double DiamHoleInLongTube    = (*m_Abrt)[0]->getFloat("CRYATTD1") * mm; //see up...
131    const double lengthTubeVoussAtt    = (*m_Abrt)[0]->getFloat("CRYATTXH") * mm;
132    const double thicknessTubeVoussAtt = (*m_Abrt)[0]->getFloat("CRYATTTH") * mm;
133    const double PosEdgeTubeVoussAtt   = (*m_Abrt)[0]->getFloat("CRYATTRX") * mm;
134    const double lengthConeVoussAtt    = (*m_Abrt)[0]->getFloat("CRYATTXS") * mm;
135    const double CryAtWiY              = (*m_Abrt)[0]->getFloat("CRYATWIY") * mm;
136    const double CryAtWXp              = (*m_Abrt)[0]->getFloat("CRYATWXP") * mm;
137    const double CryAtWXn              = (*m_Abrt)[0]->getFloat("CRYATWXN") * mm;
138    const double CryAtWBo              = (*m_Abrt)[0]->getFloat("CRYATWBO") * mm;
139    const double CryAtWTh              = (*m_Abrt)[0]->getFloat("CRYATWTH") * mm;
140    const double CryAtWZe              = (*m_Abrt)[0]->getFloat("CRYATWZE") * mm;
141    const double CryAtWRa              = (*m_Abrt)[0]->getFloat("CRYATWRA") * mm;
142    const double CryAtWYc              = (*m_Abrt)[0]->getFloat("CRYATWYC") * mm;
143    
144 
145 //!!! Fe50 (diluted iron, rho=1/2*rho_fe) ia seichas ispol'zuiu iron                                                                              //
146    const std::string voussAttMaterial = getMaterial("Iron"); 
147    const std::string voussoirMaterial = getMaterial("Aluminium");
148 //----------------------------------------------------------------------------//
149 
150    int nBTVouss = 8;                                               //
151    double zPosBTVouss[maxDim];                                                
152    for ( int i = 0; i < nBTVouss; i++ )  zPosBTVouss[i] = (*m_Abrt)[0]->getFloat("ZVOUSS",i) * mm; 
153                                                                               
154    //never used   const double zWidthBTVouss = (*m_Abrt)[0]->getFloat("VOUBLZLE") * mm;                            
155                                                                               
156 //----------------------------------------------------------------------------//
157   
158   const double tolerance = 0.5 * mm;//50
159 
160   const double phi0    = M_PI/8;
161   const double sinPhi0 = sin(phi0);
162   const double cosPhi0 = cos(phi0);
163   const double tanPhi0 = tan(phi0);
164 
165   // array of phi values (centres of phi sectors)
166   const int nPhiSectors = 8;
167   double phiVal[nPhiSectors];
168   for ( int i = 0; i < nPhiSectors; i++ )  phiVal[i] = ( 2*i + 1 ) * phi0;  
169 
170   // auxiliary arrays to store phi values, positions, rotation angles for all volumes to be replicated 
171   const int nPosMax = maxDim;
172   double phiAux[ nPhiSectors * nPosMax ], pos1Aux[ nPhiSectors * nPosMax ], pos2Aux[ nPhiSectors * nPosMax ], 
173          rotAngleAux[ nPhiSectors * nPosMax ], pos3Aux[ nPhiSectors * nPosMax];
174 
175   //------------------------------------------------------------------------------------------
176   //  Build BT Vacuum Vessel tubes
177   //------------------------------------------------------------------------------------------
178 
179 
180   // Vacuum Vessel properties 
181   const double distToCorner = sqrt(2) * rCurvBTVessel * tanPhi0; 
182   const double rIn          = radiusBTVessel - thicknessBTVessel; 
183   const double rOut         = radiusBTVessel; 
184   const double distToOut    = sqrt(2) * ( rCurvBTVessel + rOut ) * tanPhi0; 
185   const double zMaxL        = 2 * ( zMaxBTVessel - distToOut ); 
186   const double deltaR       = rMaxBTVessel - rMinBTVessel;
187   const double zMaxS        = deltaR - 2 * distToOut; 
188   const double zMaxC        = 2 * ( rCurvBTVessel + rOut ) * tanPhi0; 
189   const double rMean        = 0.5 * ( rMinBTVessel + rMaxBTVessel );  
190   
191   // Hole diameter at long tube & voussoir attachments properties 
192   const double HoleDiam         = DiamHoleOutLongTube;
193   const double HoleDiamIn       = DiamHoleInLongTube;
194   const double lengthTubeVAt    = lengthTubeVoussAtt;
195   const double thicknessTubeVAt = thicknessTubeVoussAtt;
196   const double PosEdgeTubeVAt   = PosEdgeTubeVoussAtt;
197   const double lengthConeVAt    = lengthConeVoussAtt;
198       
199   // positioning variables
200   const double xLong   = deltaR/2 - radiusBTVessel;
201   const double zShort  = zMaxBTVessel - radiusBTVessel;
202   const double xCorner = xLong  - distToCorner/2;
203   const double zCorner = zShort - distToCorner/2;
204 
205   //------------------------------------------------------------------------------------------
206   // volumes for boolean subtraction:
207 
208   // (large enough) empty box to be subtracted from sides
209   const double deltaXBox = 2 * rOut/cosPhi0;
210   const double deltaYBox = 2 * rOut;
211   const double deltaZBox = 2 * rOut * sinPhi0;
212   const double xPosBox   = rOut;  
213   
214   GeoBox* sEmptyBox = new GeoBox( deltaXBox , deltaYBox, deltaZBox );
215 
216   HepRotateY3D    rotLeftBox(-phi0), rotRightBox(phi0);
217   HepTranslateX3D trlBox(xPosBox);
218   
219   // (large enough) empty box to cut out the Rib clearances
220   GeoBox* sEmptyRibBox = new GeoBox( rOut/4, heightBTRib/2, zWidthBTRib/2 );
221 
222   HepTransform3D trlEmptyRibBox[7];
223   for (int i=0; i < nBTRibs; i++)
224   {
225   trlEmptyRibBox[i] = HepTranslateZ3D(zPosBTRib[i]) * HepTranslateX3D(rOut);
226   }
227   
228   GeoTube* sEmptyTubeVo = new GeoTube( 0., HoleDiam/2, rOut );
229 
230   HepTransform3D trlEmptyTubeVo[8];
231   for (int i=0; i < 8; i++)
232   {
233   trlEmptyTubeVo[i] = HepTranslateZ3D(zPosBTVouss[i]) * HepTranslateX3D(rOut)
234                       * HepRotateY3D(M_PI/2);
235   }
236   
237   GeoTube* sEmptyTubeCryR = new GeoTube( 0, radiusBTCryoring, rOut );
238     
239   HepTransform3D trlEmptyTubeCryR1 = HepTranslate3D( rOut * sinPhi0, rOut * cosPhi0, zPosBTCryoring ) * 
240                                      HepRotateY3D(M_PI/2) * HepRotateX3D(-M_PI/2 + phi0);
241   HepTransform3D trlEmptyTubeCryR2 = HepTranslate3D(rOut * sinPhi0, -rOut * cosPhi0, zPosBTCryoring ) *
242                                      HepRotateY3D(M_PI/2) * HepRotateX3D(M_PI/2 - phi0);
243   
244   //------------------------------------------------------------------------------------------
245   // Build long upper and lower Vacuum Vessel tubes
246   //------------------------------------------------------------------------------------------
247 
248   const GeoShape* sLongTubeIn = new GeoTube( rIn, rOut, zMaxL/2 + tolerance ); 
249 
250   const GeoShape* sLongTubeOut= new GeoTube( rIn, rOut, zMaxL/2 + tolerance );
251                                                                       
252   //------------------------------------------------------------------------------------------
253   // cut out for Voussoir ....
254   for ( int i = 0; i < nBTVouss; i++ )
255   {
256     sLongTubeIn = &( sLongTubeIn->subtract( (*sEmptyTubeVo) << trlEmptyTubeVo[i] ) );
257   } 
258   //cut out for ribs 
259   for (int i = 0; i < nBTRibs; i++)
260   {
261     sLongTubeIn = &( sLongTubeIn->subtract( (*sEmptyRibBox) << trlEmptyRibBox[i] ) ); 
262   } 
263   for (int i = 0; i < nBTRibs; i++)
264   {
265     sLongTubeOut = &( sLongTubeOut->subtract( (*sEmptyRibBox) << trlEmptyRibBox[i] ) );
266   }
267   
268   //------------------------------------------------------------------------------------------
269   // cut away left and right sides
270   HepTransform3D trlLeft  = trlBox * HepTranslateZ3D(-zMaxL/2), 
271                  trlRight = trlBox * HepTranslateZ3D( zMaxL/2);
272   
273   const GeoShape& sBevelledLongTubeIn = sLongTubeIn->subtract( (*sEmptyBox) << (trlLeft  * rotLeftBox) ).
274                                                      subtract( (*sEmptyBox) << (trlRight * rotRightBox) );
275   const GeoShape& sBevelledLongTubeOut= sLongTubeOut->subtract( (*sEmptyBox) << (trlLeft * rotLeftBox) ).
276                                                       subtract( (*sEmptyBox) << (trlRight * rotRightBox) ).
277                                                       subtract( (*sEmptyTubeCryR) << trlEmptyTubeCryR1 ).
278                                                       subtract( (*sEmptyTubeCryR) << trlEmptyTubeCryR2 );
279 
280 //inside long tube
281   GeoLogVol*  lBevelledLongTubeIn = new GeoLogVol( "BTBevelledLongTubeIn", 
282                                                  &sBevelledLongTubeIn, 
283                                                  theMaterialManager->getMaterial(vesselMaterial) );
284   GeoPhysVol* pBevelledLongTubeIn = new GeoPhysVol(lBevelledLongTubeIn);
285     
286 //outside long tube  
287   GeoLogVol*  lBevelledLongTubeOut = new GeoLogVol( "BTBevelledLongTubeOut",
288                                                  &sBevelledLongTubeOut,
289                                                  theMaterialManager->getMaterial(vesselMaterial) );
290   GeoPhysVol* pBevelledLongTubeOut = new GeoPhysVol(lBevelledLongTubeOut);
291 
292   //------------------------------------------------------------------------------------------
293   // replicate and position
294   const int nPosLongTube = 1;
295   double xPosLongTubeIn[nPosLongTube]  = { rMean - xLong };  // x positions
296   double xPosLongTubeOut[nPosLongTube] = { rMean + xLong };
297   double angleLongTube[nPosLongTube] = { M_PI };                       // rotation angles (around z)
298     
299   // fill auxiliary arrays 
300   for ( int i = 0; i < nPhiSectors * nPosLongTube; i++ )  
301   {
302      int ii = i % nPhiSectors, 
303          jj = i / nPhiSectors;
304      phiAux[i]      = phiVal[ii]; 
305      pos1Aux[i]     = xPosLongTubeIn[jj];
306      pos3Aux[i]     = xPosLongTubeOut[jj];
307      rotAngleAux[i] = angleLongTube[jj];
308   }
309 
310   GENFUNCTION fLongTube       = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosLongTube ); 
311   GENFUNCTION fTrlLongTubeIn  = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosLongTube ); 
312   GENFUNCTION fTrlLongTubeOut = ArrayFunction( pos3Aux,     pos3Aux     + nPhiSectors * nPosLongTube );
313   GENFUNCTION fRotLongTube = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosLongTube ); 
314 
315   TRANSFUNCTION XFLongTubeIn = Pow( HepRotateZ3D(1.0), fLongTube ) * 
316                                Pow( HepTranslateX3D(1.0), fTrlLongTubeIn );
317 
318   TRANSFUNCTION XFLongTubeOut= Pow( HepRotateZ3D(1.0), fLongTube ) * 
319                                Pow( HepTranslateX3D(1.0), fTrlLongTubeOut ) *
320                                Pow (HepRotateZ3D(1.0), fRotLongTube );
321                                
322   GeoSerialTransformer* sxLongTubeIn = new GeoSerialTransformer( pBevelledLongTubeIn,
323                                                                  &XFLongTubeIn,
324                                                                  nPhiSectors * nPosLongTube );
325                                                                  
326   GeoSerialTransformer* sxLongTubeOut= new GeoSerialTransformer( pBevelledLongTubeOut,
327                                                                  &XFLongTubeOut,
328                                                                  nPhiSectors * nPosLongTube );
329   container->add(sxLongTubeIn);
330   container->add(sxLongTubeOut);
331 
332 //--------------VoussoirAttWing---------
333 
334   GeoShape* sVoussoirAttWingBox = new GeoBox( CryAtWiY, CryAtWZe/2 + CryAtWTh, ( CryAtWXp - CryAtWXn
335                                               - CryAtWBo - tolerance )/2 );
336   GeoTrd* sVoussoirAttWingTrap1 = new GeoTrd( CryAtWiY, CryAtWiY - CryAtWBo, CryAtWZe/2 + CryAtWTh,
337                                               CryAtWZe/2 + CryAtWTh, CryAtWBo/2 );
338   GeoTrd* sVoussoirAttWingTrap2 = new GeoTrd( CryAtWiY -( 2*radiusBTVessel + CryAtWXn )
339                                               *sin(phi0), CryAtWiY, CryAtWZe/2 + CryAtWTh, CryAtWZe/2
340                                               + CryAtWTh, radiusBTVessel + CryAtWXn/2 - tolerance );
341   GeoBox* sEmptyVoussoirAttWingBox1 = new GeoBox( 2*radiusBTVessel, CryAtWZe/2, 3*radiusBTVessel );
342   GeoBox* sEmptyVoussoirAttWingBox2 = new GeoBox( CryAtWYc, CryAtWZe, radiusBTVessel );
343   GeoTubs* sEmptyVoussoirAttWingTube1 = new GeoTubs( CryAtWRa, 2*CryAtWRa, CryAtWZe, M_PI + phi0,
344                                                    6 * phi0 );
345   GeoTube* sEmptyVoussoirAttWingTube2 = new GeoTube( 0, radiusBTVessel + tolerance, CryAtWZe );
346 
347   const GeoShape& sVoussoirAttWing = sVoussoirAttWingBox->add( (*sVoussoirAttWingTrap1) << 
348                                                           HepTranslateZ3D( ( CryAtWXp - 
349                                                           CryAtWXn )/2 ) ).
350                                                           add( (*sVoussoirAttWingTrap2) <<
351                                                           HepTranslateZ3D( -( ( CryAtWXp - 
352                                                           CryAtWBo )/2
353                                                           + radiusBTVessel ) ) ).
354                                                           subtract( (*sEmptyVoussoirAttWingBox1) << 
355                                                           HepTranslateZ3D( ( CryAtWBo - CryAtWXp - 
356                                                           CryAtWXn )/2 ) ).
357                                                           subtract( (*sEmptyVoussoirAttWingBox2) << 
358                                                           HepTranslateZ3D( ( CryAtWBo - CryAtWXp -
359                                                           CryAtWXn )/2 
360                                                           - radiusBTVessel )  ).
361                                                           subtract( (*sEmptyVoussoirAttWingTube1) <<
362                                                           HepTranslateZ3D( ( CryAtWBo - CryAtWXp -
363                                                           CryAtWXn )/2 ) * HepRotateX3D(M_PI/2) ).
364                                                           subtract( (*sEmptyVoussoirAttWingTube2) <<
365                                                           HepTranslateZ3D( ( CryAtWBo - CryAtWXp -
366                                                           CryAtWXn )/2 ) * HepRotateX3D(M_PI/2) );
367 
368   GeoLogVol* lVoussoirAttWing = new GeoLogVol( "BTVoussoirAttWing", &sVoussoirAttWing,
369                                                theMaterialManager->getMaterial(voussoirMaterial) );
370   GeoPhysVol* pVoussoirAttWing = new GeoPhysVol(lVoussoirAttWing);
371 
372   // replicate and position
373   for (int i = 0; i < nPhiSectors * nBTVouss; i++ )
374   {
375     int ii = i % nPhiSectors,
376         jj = i / nPhiSectors;
377     phiAux[i] = phiVal[ii];
378     pos1Aux[i] = zPosBTVouss[jj];
379   }
380 
381   GENFUNCTION fVoussoirAttWing    = ArrayFunction( phiAux, phiAux + nPhiSectors * nBTVouss );
382   GENFUNCTION fTrlVoussoirAttWing = ArrayFunction( pos1Aux, pos1Aux + nPhiSectors * nBTVouss );
383 
384   TRANSFUNCTION XFVoussoirAttWing = Pow( HepRotateZ3D(1.0), fVoussoirAttWing ) *
385                               Pow( HepTranslateZ3D(1.0), fTrlVoussoirAttWing ) *
386                               HepTranslateX3D( rMean - (CryAtWBo - CryAtWXp - CryAtWXn)/2 - xLong )
387                               * HepRotateX3D(M_PI/2) * HepRotateY3D(M_PI/2);
388   GeoSerialTransformer* sxVoussoirAttWing = new GeoSerialTransformer( pVoussoirAttWing,
389                                                                 &XFVoussoirAttWing,
390                                                                 nPhiSectors * nBTVouss );
391   container->add(sxVoussoirAttWing);
392 
393   //------------------------------------------------------------------------------------------
394   // Build short left and right Vacuum Vessel tubes
395   //------------------------------------------------------------------------------------------
396 
397   GeoTube* sShortTube = new GeoTube( rIn, rOut, zMaxS/2 - tolerance );
398 
399   //------------------------------------------------------------------------------------------
400   // cut away left and right sides
401   trlLeft  = trlBox * HepTranslateZ3D(-zMaxS/2); 
402   trlRight = trlBox * HepTranslateZ3D( zMaxS/2);
403   const GeoShape& sBevelledShortTube = sShortTube->subtract( (*sEmptyBox) << (trlLeft  * rotLeftBox) ).
404                                                    subtract( (*sEmptyBox) << (trlRight * rotRightBox) );
405   GeoLogVol*  lBevelledShortTube = new GeoLogVol( "BTBevelledShortTube", 
406                                                   &sBevelledShortTube, 
407                                                   theMaterialManager->getMaterial(vesselMaterial) );
408   GeoPhysVol* pBevelledShortTube = new GeoPhysVol(lBevelledShortTube);
409 
410   //------------------------------------------------------------------------------------------
411   // replicate and position
412   const int nPosShortTube = 2;
413   double zPosShortTube[nPosShortTube]  = { -zShort, zShort };  // z positions
414   double angleShortTube[nPosShortTube] = { -M_PI/2, M_PI/2 };  // rotation angles (around y)
415     
416   // fill auxiliary arrays 
417   for ( int i = 0; i < nPhiSectors * nPosShortTube; i++ )  
418   {
419      int ii = i % nPhiSectors, 
420          jj = i / nPhiSectors;
421      phiAux[i]      = phiVal[ii]; 
422      pos1Aux[i]     = zPosShortTube[jj];
423      rotAngleAux[i] = angleShortTube[jj];
424   }
425 
426   GENFUNCTION fShortTube    = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosShortTube ); 
427   GENFUNCTION fTrlShortTube = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosShortTube ); 
428   GENFUNCTION fRotShortTube = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosShortTube ); 
429 
430   TRANSFUNCTION XFShortTube = Pow( HepRotateZ3D(1.0), fShortTube ) * 
431                               Pow( HepTranslateZ3D(1.0), fTrlShortTube ) * HepTranslateX3D(rMean) *
432                               Pow( HepRotateY3D(1.0), fRotShortTube );
433   GeoSerialTransformer* sxShortTube = new GeoSerialTransformer( pBevelledShortTube, 
434                                                                 &XFShortTube, 
435                                                                 nPhiSectors * nPosShortTube );
436   container->add(sxShortTube);
437 
438 
439   //------------------------------------------------------------------------------------------
440   // Build corner Vacuum Vessel tubes
441   //------------------------------------------------------------------------------------------
442 
443   GeoTube* sCornerTube = new GeoTube( rIn, rOut - tolerance, zMaxC/2 );
444 
445   //------------------------------------------------------------------------------------------
446   // cut away left and right sides
447   trlLeft  = trlBox * HepTranslateZ3D(-zMaxC/2); 
448   trlRight = trlBox * HepTranslateZ3D( zMaxC/2);
449   const GeoShape& sBevelledCornerTube = sCornerTube->subtract( (*sEmptyBox) << (trlLeft  * rotLeftBox) ).
450                                                      subtract( (*sEmptyBox) << (trlRight * rotRightBox) );
451   GeoLogVol*  lBevelledCornerTube = new GeoLogVol( "BTBevelledCornerTube", 
452                                                    &sBevelledCornerTube,
453                                                    theMaterialManager->getMaterial(vesselMaterial) );
454   GeoPhysVol* pBevelledCornerTube = new GeoPhysVol(lBevelledCornerTube);
455 
456   //------------------------------------------------------------------------------------------
457   // replicate and position
458   const int nPosCornerTube = 4;
459   double zPosCornerTube[nPosCornerTube]  = { zCorner, zCorner, -zCorner, -zCorner };  // z positions
460   double xPosCornerTube[nPosCornerTube]  = { rMean - xCorner, rMean + xCorner,  
461                                              rMean + xCorner, rMean - xCorner };      // x positions
462   double angleCornerTube[nPosCornerTube] = { M_PI/4, 3*M_PI/4, 5*M_PI/4, 7*M_PI/4 };  // rotation angles (around y)
463     
464   // fill auxiliary arrays 
465   for ( int i = 0; i < nPhiSectors * nPosCornerTube; i++ )  
466   {
467      int ii = i % nPhiSectors, 
468          jj = i / nPhiSectors;
469      phiAux[i]      = phiVal[ii]; 
470      pos1Aux[i]     = zPosCornerTube[jj];
471      pos2Aux[i]     = xPosCornerTube[jj];
472      rotAngleAux[i] = angleCornerTube[jj];
473   }
474 
475   GENFUNCTION fCornerTube     = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosCornerTube ); 
476   GENFUNCTION fTrl1CornerTube = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosCornerTube ); 
477   GENFUNCTION fTrl2CornerTube = ArrayFunction( pos2Aux,     pos2Aux     + nPhiSectors * nPosCornerTube ); 
478   GENFUNCTION fRotCornerTube  = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosCornerTube ); 
479 
480   TRANSFUNCTION XFCornerTube = Pow( HepRotateZ3D(1.0), fCornerTube ) * 
481                                Pow( HepTranslateZ3D(1.0), fTrl1CornerTube ) *
482                                Pow( HepTranslateX3D(1.0), fTrl2CornerTube ) *
483                                Pow( HepRotateY3D(1.0), fRotCornerTube );
484   GeoSerialTransformer* sxCornerTube = new GeoSerialTransformer( pBevelledCornerTube,  
485                                                                  &XFCornerTube, 
486                                                                  nPhiSectors * nPosCornerTube );   
487   container->add(sxCornerTube);
488   
489   
490   //------------------------------------------------------------------------------------------
491   //  Build BT Vacuum Vessel ribs
492   //  - account for correct concave profile at sides touching the Vacuum Vessel tubes
493   //------------------------------------------------------------------------------------------
494    
495   // rib property
496   const double lengthBTRib = deltaR - 4 * rOut;
497   
498   //------------------------------------------------------------------------------------------
499   // volumes for boolean subtraction:
500 
501   // (large enough) empty tube to cut out the concave profile from the sides
502   double alpha   = asin( 0.5 * heightBTRib/rOut );
503   double sagitta = rOut - 0.5 * heightBTRib/tan(alpha);
504 
505   GeoShape* sEmptyTube = new GeoTube( 0, rOut, 2*zWidthBTRib - tolerance );
506   
507   // (large enough) empty box to cut out the inner Rib volume
508   GeoBox* sEmptyInnerBox = new GeoBox( lengthBTRib/2 + rOut/2 , 
509                                        heightBTRib/2 - thicknessBTVessel, 
510                                        zWidthBTRib/2 - thicknessBTVessel );
511 
512   //------------------------------------------------------------------------------------------
513   // rib envelope
514   GeoBox* sRibOuterBox = new GeoBox( lengthBTRib/2 + sagitta, heightBTRib/2, zWidthBTRib/2 ); 
515   
516   // cut away concave profile on left and right sides, and cut out inner volume
517   const GeoShape& sRibEnvelope = sRibOuterBox->subtract( (*sEmptyTube) << HepTranslateX3D( -lengthBTRib/2 - rOut ) ).
518                                                subtract( (*sEmptyTube) << HepTranslateX3D(  lengthBTRib/2 + rOut ) ).
519                                                subtract( *sEmptyInnerBox );
520 
521   GeoLogVol*  lRibEnvelope = new GeoLogVol( "BTRibEnvelope", 
522                                             &sRibEnvelope, 
523                                             theMaterialManager->getMaterial(vesselMaterial) );
524   GeoPhysVol* pRibEnvelope = new GeoPhysVol(lRibEnvelope);
525 
526   //------------------------------------------------------------------------------------------
527   // replicate and position
528   for ( int i = 0; i < nPhiSectors * nBTRibs; i++ )  // fill auxiliary arrays 
529   {
530      int ii = i % nPhiSectors, 
531          jj = i / nPhiSectors;
532      phiAux[i]  = phiVal[ii]; 
533      pos1Aux[i] = zPosBTRib[jj];
534   }
535 
536   GENFUNCTION fRibEnvelope    = ArrayFunction( phiAux,  phiAux  + nPhiSectors * nBTRibs ); 
537   GENFUNCTION fTrlRibEnvelope = ArrayFunction( pos1Aux, pos1Aux + nPhiSectors * nBTRibs ); 
538 
539   TRANSFUNCTION XFRibEnvelope = Pow( HepRotateZ3D(1.0), fRibEnvelope ) * 
540                                 Pow( HepTranslateZ3D(1.0), fTrlRibEnvelope ) 
541                                 * HepTranslateX3D(rMean);
542   GeoSerialTransformer* sxRibEnvelope = new GeoSerialTransformer( pRibEnvelope, 
543                                                                   &XFRibEnvelope, 
544                                                                   nPhiSectors * nBTRibs );
545   container->add(sxRibEnvelope);
546 
547 //Build BT Vacuum Vessel wings ribs
548 // DHW - June 2008: temporarily take 25 mm off of CryRiWTh
549   double dhw = 25.;
550   GeoShape* sTrapWingRib   = new GeoTrd( CryRiWYn/2, CryRiWYp/2, zWidthBTRib/2 + CryRiWTh - dhw,
551                                         zWidthBTRib/2 + CryRiWTh - dhw, (CryRiWXp - CryRiWXn)/2 );
552   GeoBox* sBoxWingRib      = new GeoBox( CryRiWYn/2, zWidthBTRib/2 + CryRiWTh - dhw,
553                                         (rOut + CryRiWXn)/2 );
554   GeoBox* sEmptyBoxWingRib = new GeoBox( 2 * rOut, zWidthBTRib/2 + tolerance, 2 * rOut );
555   GeoTube* sEmptyTubeWingRib= new GeoTube( 0, rOut + tolerance, zWidthBTRib);
556   
557   const GeoShape& sWingRib = sTrapWingRib->add( (*sBoxWingRib) << HepTranslateZ3D(
558                                            -(rOut + CryRiWXp)/2) ).
559                                            subtract( (*sEmptyBoxWingRib) << HepTranslateZ3D(
560                                            -(CryRiWXp + CryRiWXn)/2) ).
561                                            subtract( (*sEmptyTubeWingRib) << HepTranslateZ3D(
562                                            -(CryRiWXp + CryRiWXn + 2*rOut)/2) * 
563                                            HepRotateX3D(M_PI/2) );
564   
565   GeoLogVol* lWingRib = new GeoLogVol( "BTWingRib", &sWingRib,
566                                         theMaterialManager->getMaterial(voussoirMaterial) );
567   GeoPhysVol* pWingRib = new GeoPhysVol(lWingRib);
568                                         
569   // replicate and position
570   for (int i = 0; i < nPhiSectors * nBTRibs; i++ )
571   {
572     int ii = i % nPhiSectors,
573         jj = i / nPhiSectors;
574     phiAux[i] = phiVal[ii];
575     pos1Aux[i] = zPosBTRib[jj];
576   }
577   
578   GENFUNCTION fWingRib    = ArrayFunction( phiAux, phiAux + nPhiSectors * nBTRibs );
579   GENFUNCTION fTrlWingRib = ArrayFunction( pos1Aux, pos1Aux + nPhiSectors * nBTRibs );
580   
581   TRANSFUNCTION XFWingRibIn = Pow( HepRotateZ3D(1.0), fWingRib ) * 
582                               Pow( HepTranslateZ3D(1.0), fTrlWingRib ) *
583                               HepTranslateX3D( rMinBTVessel + 2*rOut + (CryRiWXp + CryRiWXn)/2)
584                               * HepRotateX3D(M_PI/2) * HepRotateY3D(M_PI/2);
585   GeoSerialTransformer* sxWingRibIn = new GeoSerialTransformer( pWingRib,
586                                                                 &XFWingRibIn,
587                                                                 nPhiSectors * nBTRibs );
588   TRANSFUNCTION XFWingRibOut = Pow( HepRotateZ3D(1.0), fWingRib )*
589                                Pow( HepTranslateZ3D(1.), fTrlWingRib ) * 
590                                HepTranslateX3D( rMaxBTVessel - 2*rOut - (CryRiWXp + CryRiWXn)/2)
591                                * HepRotateX3D(M_PI/2) * HepRotateY3D(-M_PI/2);
592   GeoSerialTransformer* sxWingRibOut = new GeoSerialTransformer( pWingRib,
593                                                                  &XFWingRibOut,
594                                                                  nPhiSectors * nBTRibs );
595   container->add(sxWingRibIn);                                                        
596   container->add(sxWingRibOut);
597   //------------------------------------------------------------------------------------------
598   //  Build BT Voussoir Attachments
599   //  - trapezoidal shapes, embracing the long tubes; material: diluted Iron (approximation)
600   //  - account for correct concave profile at sides touching the Vacuum Vessel tubes
601   //------------------------------------------------------------------------------------------
602   
603   GeoTube* sVoussAttTube  = new GeoTube( HoleDiam/2 - thicknessTubeVAt, HoleDiam/2,
604                                         (lengthTubeVAt - lengthConeVAt)/2 );
605 //never used   GeoCons* sVoussAttCone1 = new GeoCons( HoleDiam/2 - thicknessTubeVAt, HoleDiamIn/2 -thicknessTubeVAt,
606 //never used                                     HoleDiam/2, HoleDiamIn/2,
607 //never used                                     lengthConeVAt/2, 0., M_PI );
608   GeoCons* sVoussAttCone2 = new GeoCons( HoleDiam/2 - thicknessTubeVAt, HoleDiamIn/2 - thicknessTubeVAt,
609                                          HoleDiam/2, HoleDiamIn/2,
610                                          lengthConeVAt/2, 0, 2*M_PI);
611 
612   //------------------------------------------------------------------------------------------
613   // cut out concave profile on bottom side 
614   // - note that empty tube axes don't coincide with trapezoid ones
615   double xPosVoussAtt = PosEdgeTubeVAt - lengthTubeVAt/2 - lengthConeVAt/2;
616   
617   HepTransform3D xfEmptyTube = HepTranslateZ3D( -( xPosVoussAtt + xLong - rMean ) )
618                                                 * HepRotateX3D( M_PI/2 );
619   const GeoShape& sVoussAtt  = sVoussAttTube->subtract( (*sEmptyTube) << xfEmptyTube ).
620 //                                            add( (*sVoussAttCone1) << HepTranslateZ3D( lengthTubeVAt/2 ) ).
621                                               add( (*sVoussAttCone2) << HepTranslateZ3D( lengthTubeVAt/2 ) );
622                                             
623   GeoLogVol*  lVoussAtt = new GeoLogVol( "BTVoussoirAttachment", 
624                                          &sVoussAtt, 
625                                          theMaterialManager->getMaterial(voussAttMaterial) );
626   GeoPhysVol* pVoussAtt = new GeoPhysVol(lVoussAtt);
627 
628   //------------------------------------------------------------------------------------------
629   // replicate and position
630  
631 //  double xPosVoussAtt = PosEdgeTubeVAt - lengthTubeVAt/2 - lengthConeVAt/2;
632   
633   for ( int i = 0; i < nPhiSectors * nBTVouss; i++ )  // fill auxiliary arrays 
634   {
635      int ii = i % nPhiSectors, 
636          jj = i / nPhiSectors;
637      phiAux[i]  = phiVal[ii]; 
638      pos1Aux[i] = zPosBTVouss[jj];
639   }
640 
641   GENFUNCTION fVoussAtt     = ArrayFunction( phiAux,  phiAux  + nPhiSectors * nBTVouss ); 
642   GENFUNCTION fTrlVoussAtt1 = ArrayFunction( pos1Aux, pos1Aux + nPhiSectors * nBTVouss ); 
643 
644   TRANSFUNCTION XFVoussAtt = Pow( HepRotateZ3D(1.0), fVoussAtt ) * 
645                              Pow( HepTranslateZ3D(1.0), fTrlVoussAtt1 ) *
646                              HepTranslateX3D(xPosVoussAtt) *
647                              HepRotateX3D( -M_PI/2 ) * HepRotateY3D(M_PI/2);
648   GeoSerialTransformer* sxVoussAtt = new GeoSerialTransformer( pVoussAtt, &XFVoussAtt, nPhiSectors * nBTVouss );
649   container->add(sxVoussAtt);
650 
651 }
652  
653 void BarrelToroidBuilderRDB::buildBTColdMass( GeoPhysVol* container ) 
654 {
655   ////////////////////////////////////////////////////////////////////////////////////////////
656   //  Builds BT Cold Mass (coils, coil casings, and ribs)                                   //
657   ////////////////////////////////////////////////////////////////////////////////////////////
658 
659   const StoredMaterialManager*  theMaterialManager;
660   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
661   {
662     return;
663   } 
664 
665   const int maxDim = 7; 
666 
667 
668 //---------- BTCB -- BT Conductor and Coil Casing Parameters -----------------//
669 //                                                                            //
670                                                                               //
671    const double widthBTColdMass  = (*m_Abrt)[0]->getFloat("COMARTHI") * mm;                                 //btcb->xwidth
672    const double heightBTColdMass = (*m_Abrt)[0]->getFloat("COMAYTHI") * mm;                                 //btcb->ywidth
673                                                                               // 
674    const std::string conductorMaterial = getMaterial("Aluminium"); 
675 //                                                                            //
676 //----------------------------------------------------------------------------//
677   
678 //---------- BTVV -- BT Vacuum Vessel Parameters -----------------------------//
679 //                                                                            //
680                                                                               //
681    const double rMinBTVessel      = (*m_Abrt)[0]->getFloat("COMARMIN") * mm;  //CryoRmin/101
682    const double rMaxBTVessel      = (*m_Abrt)[0]->getFloat("COMARMAX") * mm;  //CryiRmax/102
683    const double zMaxBTVessel      = (*m_Abrt)[0]->getFloat("COMAZMAX") * mm;  //CryoZmax/103
684    const double rCurvBTVessel     = (*m_Abrt)[0]->getFloat("COMARCUI") * mm;  //CryoRcurv/104
685 //   const double radiusBTVessel    = (*m_Abrt)[0]->getFloat("CRYORADI") * mm;  //  
686 //----------------------------------------------------------------------------//
687 
688 //---------- BTRC -- BT Rib and Rib Cold Mass Parameters ---------------------//
689 //                                                                            //
690                                                                               //
691    int nBTRibs = 7;                                                           //btrc->nrib
692    double zPosBTRib[maxDim];                                                  //
693    for ( int i = 0; i < nBTRibs; i++ )  zPosBTRib [i] = (*m_Abrt)[0]->getFloat("ZRIB",i) * mm; //
694                                                                               // 
695    const double heightBTRibColdMass = (*m_Abrt)[0]->getFloat("COMARIBY") * mm;//ColdMassRibY/173
696    const double zWidthBTRibColdMass = (*m_Abrt)[0]->getFloat("COMARIBZ") * mm;// ColdMassRibZ/174
697                                                                               //
698 //----------------------------------------------------------------------------//
699   
700 
701   const double phi0    = M_PI/8;
702 //never used   const double sinPhi0 = sin(phi0);
703 //never used   const double cosPhi0 = cos(phi0);
704   const double tanPhi0 = tan(phi0);
705 
706   // array of phi values (centres of phi sectors)
707   const int nPhiSectors = 8;
708   double phiVal[nPhiSectors];
709   for ( int i = 0; i < nPhiSectors; i++ )  phiVal[i] = ( 2*i + 1 ) * phi0;  
710 
711   // auxiliary arrays to store phi values, positions, rotation angles for all volumes to be replicated 
712   const int nPosMax = maxDim;
713   double phiAux[ nPhiSectors * nPosMax ], pos1Aux[ nPhiSectors * nPosMax ], pos2Aux[ nPhiSectors * nPosMax ], 
714          rotAngleAux[ nPhiSectors * nPosMax ];
715 
716 
717   // Cold Mass properties 
718   const double deltaL       = widthBTColdMass * tanPhi0; 
719   const double distToCorner = sqrt(2) * rCurvBTVessel * tanPhi0; 
720   const double rMean        = 0.5 * ( rMinBTVessel + rMaxBTVessel );  
721   const double deltaR       = rMaxBTVessel - rMinBTVessel;
722   const double lMeanLong    = 2 * ( zMaxBTVessel - widthBTColdMass/2 - distToCorner ); 
723   const double lMeanShort   = deltaR - 2 * ( widthBTColdMass/2 + distToCorner ); 
724   const double lMeanCorner  = 2 * rCurvBTVessel * tanPhi0; 
725 
726   // positioning variables
727   const double xLong   = deltaR/2 - widthBTColdMass/2;
728   const double zShort  = zMaxBTVessel - widthBTColdMass/2;
729   const double xCorner = xLong  - distToCorner/2;
730   const double zCorner = zShort - distToCorner/2;
731 
732 
733   //------------------------------------------------------------------------------------------
734   // Build BT Coil+Coil Casing, using trapezoidal shapes everywhere as an approximation. 
735   // Note that GeoTrd axes definition does not coincide with the coordinates used here.   
736   //  --> Shapes need to be rotated first, then put into right orientation.
737   //------------------------------------------------------------------------------------------
738 
739   //------------------------------------------------------------------------------------------
740   // Build long upper and lower Conductor Box segments 
741   //------------------------------------------------------------------------------------------
742 
743   GeoTrd* sLongSegment = new GeoTrd( lMeanLong/2 - deltaL/2, lMeanLong/2 + deltaL/2,
744                                      heightBTColdMass/2, heightBTColdMass/2, widthBTColdMass/2 );
745   
746   GeoLogVol*  lLongSegment = new GeoLogVol( "BTColdLongSegment", 
747                                             sLongSegment, 
748                                             theMaterialManager->getMaterial(conductorMaterial) );
749   GeoPhysVol* pLongSegment = new GeoPhysVol(lLongSegment);
750       
751   //------------------------------------------------------------------------------------------
752   // replicate and position
753   const int nPosLongSegment = 2;
754   double xPosLongSegment[nPosLongSegment]  = { rMean - xLong, rMean + xLong };  // x positions
755   double angleLongSegment[nPosLongSegment] = { -M_PI/2, M_PI/2 };               // rotation angles (around y)
756     
757   // fill auxiliary arrays 
758   for ( int i = 0; i < nPhiSectors * nPosLongSegment; i++ )  
759   {
760      int ii = i % nPhiSectors, 
761          jj = i / nPhiSectors;
762      phiAux[i]      = phiVal[ii]; 
763      pos1Aux[i]     = xPosLongSegment[jj];
764      rotAngleAux[i] = angleLongSegment[jj];
765   }
766 
767   GENFUNCTION fLongSegment    = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosLongSegment ); 
768   GENFUNCTION fTrlLongSegment = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosLongSegment ); 
769   GENFUNCTION fRotLongSegment = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosLongSegment ); 
770 
771   TRANSFUNCTION XFLongSegment = Pow( HepRotateZ3D(1.0), fLongSegment ) * 
772                                 Pow( HepTranslateX3D(1.0), fTrlLongSegment ) * 
773                                 Pow( HepRotateY3D(1.0), fRotLongSegment );
774   GeoSerialTransformer* sxLongSegment = new GeoSerialTransformer( pLongSegment, 
775                                                                   &XFLongSegment, 
776                                                                   nPhiSectors * nPosLongSegment );
777 
778   container->add(sxLongSegment);
779 
780 
781   //------------------------------------------------------------------------------------------
782   // Build short left and right Conductor Box segments
783   //------------------------------------------------------------------------------------------
784 
785   GeoTrd* sShortSegment = new GeoTrd( lMeanShort/2 - deltaL/2, lMeanShort/2 + deltaL/2,
786                                       heightBTColdMass/2, heightBTColdMass/2, widthBTColdMass/2 );
787   
788   GeoLogVol*  lShortSegment = new GeoLogVol( "BTColdShortSegment", 
789                                              sShortSegment, 
790                                              theMaterialManager->getMaterial(conductorMaterial) );
791   GeoPhysVol* pShortSegment = new GeoPhysVol(lShortSegment);
792       
793   //------------------------------------------------------------------------------------------
794   // replicate and position
795   const int nPosShortSegment = 2;
796   double zPosShortSegment[nPosShortSegment]  = { -zShort, zShort };  // z positions
797   double angleShortSegment[nPosShortSegment] = { M_PI, 0 };          // rotation angles (around y)
798     
799   // fill auxiliary arrays 
800   for ( int i = 0; i < nPhiSectors * nPosShortSegment; i++ )  
801   {
802      int ii = i % nPhiSectors, 
803          jj = i / nPhiSectors;
804      phiAux[i]      = phiVal[ii]; 
805      pos1Aux[i]     = zPosShortSegment[jj];
806      rotAngleAux[i] = angleShortSegment[jj];
807   }
808 
809   GENFUNCTION fShortSegment    = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosShortSegment ); 
810   GENFUNCTION fTrlShortSegment = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosShortSegment ); 
811   GENFUNCTION fRotShortSegment = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosShortSegment ); 
812 
813   TRANSFUNCTION XFShortSegment = Pow( HepRotateZ3D(1.0), fShortSegment ) * 
814                                  Pow( HepTranslateZ3D(1.0), fTrlShortSegment ) * HepTranslateX3D(rMean) *
815                                  Pow( HepRotateY3D(1.0), fRotShortSegment );
816   GeoSerialTransformer* sxShortSegment = new GeoSerialTransformer( pShortSegment, 
817                                                                    &XFShortSegment, 
818                                                                    nPhiSectors * nPosShortSegment );
819   container->add(sxShortSegment);
820 
821 
822   //------------------------------------------------------------------------------------------
823   // Build corner Conductor Box segments
824   //------------------------------------------------------------------------------------------
825 
826   GeoTrd* sCornerSegment = new GeoTrd( lMeanCorner/2 - deltaL/2, lMeanCorner/2 + deltaL/2,
827                                        heightBTColdMass/2, heightBTColdMass/2, widthBTColdMass/2 );
828 
829   GeoLogVol*  lCornerSegment = new GeoLogVol( "BTColdCornerSegment", 
830                                               sCornerSegment, 
831                                               theMaterialManager->getMaterial(conductorMaterial) );
832   GeoPhysVol* pCornerSegment = new GeoPhysVol(lCornerSegment);
833       
834   //------------------------------------------------------------------------------------------
835   // replicate and position
836   const int nPosCornerSegment = 4;
837   double zPosCornerSegment[nPosCornerSegment]  = { -zCorner, zCorner, zCorner, -zCorner }; // z positions
838   double xPosCornerSegment[nPosCornerSegment]  = { rMean + xCorner, rMean + xCorner,  
839                                                    rMean - xCorner, rMean - xCorner };     // x positions
840   double angleCornerSegment[nPosCornerSegment] = { 3*M_PI/4, M_PI/4, -M_PI/4, -3*M_PI/4 }; // rotation angles (around y)
841     
842   // fill auxiliary arrays 
843   for ( int i = 0; i < nPhiSectors * nPosCornerSegment; i++ )  
844   {
845      int ii = i % nPhiSectors, 
846          jj = i / nPhiSectors;
847      phiAux[i]      = phiVal[ii]; 
848      pos1Aux[i]     = zPosCornerSegment[jj];
849      pos2Aux[i]     = xPosCornerSegment[jj];
850      rotAngleAux[i] = angleCornerSegment[jj];
851   }
852 
853   GENFUNCTION fCornerSegment     = ArrayFunction( phiAux,      phiAux      + nPhiSectors * nPosCornerSegment ); 
854   GENFUNCTION fTrl1CornerSegment = ArrayFunction( pos1Aux,     pos1Aux     + nPhiSectors * nPosCornerSegment ); 
855   GENFUNCTION fTrl2CornerSegment = ArrayFunction( pos2Aux,     pos2Aux     + nPhiSectors * nPosCornerSegment ); 
856   GENFUNCTION fRotCornerSegment  = ArrayFunction( rotAngleAux, rotAngleAux + nPhiSectors * nPosCornerSegment ); 
857 
858   TRANSFUNCTION XFCornerSegment = Pow( HepRotateZ3D(1.0), fCornerSegment ) * 
859                                   Pow( HepTranslateZ3D(1.0), fTrl1CornerSegment ) *
860                                   Pow( HepTranslateX3D(1.0), fTrl2CornerSegment ) *
861                                   Pow( HepRotateY3D(1.0), fRotCornerSegment );
862   GeoSerialTransformer* sxCornerSegment = new GeoSerialTransformer( pCornerSegment, 
863                                                                     &XFCornerSegment, 
864                                                                     nPhiSectors * nPosCornerSegment );
865   container->add(sxCornerSegment);
866   
867   
868   //------------------------------------------------------------------------------------------
869   //  Build BT Rib Cold Mass
870   // - Cold Mass length defined by Vacuum Vessel and Coil Casing dimensions 
871   //------------------------------------------------------------------------------------------
872    
873   const double lengthBTRib = deltaR - 2 * widthBTColdMass;
874 
875   GeoBox*     sRibColdMass = new GeoBox( lengthBTRib/2, heightBTRibColdMass/2, zWidthBTRibColdMass ); 
876   
877   GeoLogVol*  lRibColdMass = new GeoLogVol( "BTColdRib", 
878                                             sRibColdMass, 
879                                             theMaterialManager->getMaterial(conductorMaterial) );
880   GeoPhysVol* pRibColdMass = new GeoPhysVol(lRibColdMass);
881 
882   //------------------------------------------------------------------------------------------
883   // replicate and position
884   for ( int i = 0; i < nPhiSectors * nBTRibs; i++ )  // fill auxiliary arrays 
885   {
886      int ii = i % nPhiSectors, 
887          jj = i / nPhiSectors;
888      phiAux[i]  = phiVal[ii]; 
889      pos1Aux[i] = zPosBTRib[jj];
890   }
891 
892   GENFUNCTION fRibColdMass    = ArrayFunction( phiAux,  phiAux  + nPhiSectors * nBTRibs ); 
893   GENFUNCTION fTrlRibColdMass = ArrayFunction( pos1Aux, pos1Aux + nPhiSectors * nBTRibs ); 
894 
895   TRANSFUNCTION XFRibColdMass = Pow( HepRotateZ3D(1.0), fRibColdMass ) * 
896                                 Pow( HepTranslateZ3D(1.0), fTrlRibColdMass ) * HepTranslateX3D(rMean);
897   GeoSerialTransformer* sxRibColdMass = new GeoSerialTransformer( pRibColdMass, &XFRibColdMass, nPhiSectors * nBTRibs );
898   container->add(sxRibColdMass);
899 
900 }
901 
902   void BarrelToroidBuilderRDB::buildBTVoussoirs( GeoPhysVol* container ) 
903 {
904   ////////////////////////////////////////////////////////////////////////////////////////////
905   //  Builds BT Voussoirs (warm support structure)                                          //
906   ////////////////////////////////////////////////////////////////////////////////////////////
907 
908   const StoredMaterialManager*  theMaterialManager;
909   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
910   {
911     return;
912   } 
913 
914   const int maxDim = 8; 
915 
916 
917 //---------- BTVV -- BT Vacuum Vessel Parameters -----------------------------//
918 //                                                                            //
919                                                                               //
920 //never used    const double rMinBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMIN") * mm;  //CryoRmin/101
921 //never used    const double rMaxBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMAX") * mm;  //CryiRmax/102
922 //                                                                            //
923 //----------------------------------------------------------------------------//
924 
925 //---------- BTVS -- BT Voussoir Parameters ----------------------------------//
926 //                                                                            //
927                                                                               //
928    int nBTVouss = 8;                                               //btvs->nvouss
929    double zPosBTVouss[maxDim];                                                //
930    for ( int i = 0; i < nBTVouss; i++ )  zPosBTVouss[i] = (*m_Abrt)[0]->getFloat("ZVOUSS", i) * mm; //
931                                                                               //
932    const double tRMinBTVouss                 = (*m_Abrt)[0]->getFloat("VOUSMBRA") * mm;           //btvs->trmin est' srednee snach etih dvuh per. ->
933 //never used    const double tRBTVouss                    = (*m_Abrt)[0]->getFloat("VOUSSRAD") * mm; 
934    const double tWidthBTVouss                = (*m_Abrt)[0]->getFloat("VOUSBLXH") * mm;//
935    
936    const std::string voussoirMaterial = getMaterial("Aluminium"); 
937                                                                               //
938 //----------------------------------------------------------------------------//
939 
940    const double zWidth               = (*m_Abrt)[0]->getFloat("VOUBLZLE") * mm;
941    const double zWidthStepVouss       = (*m_Abrt)[0]->getFloat("VOUBLZLS") * mm;
942    const double alphaStep                     = 15.0 * deg;
943    const double xLength       = (*m_Abrt)[0]->getFloat("VOUSBLYW") * mm;
944    const double thicknessStepVouss     = (*m_Abrt)[0]->getFloat("VOURECSL") * mm;
945 //never used    const double xPosFoot                 = (*m_Feet)[0]->getFloat("STDFOOXP") * mm;
946 //never used    const double yPosFoot                 = (*m_Feet)[0]->getFloat("STDFOOYP") * mm;
947    const double distInnerRib    = 25.0 + (*m_Abrt)[0]->getFloat("VOUCRCYR") * mm;
948    const double distHoleCentre  = (*m_Abrt)[0]->getFloat("VOUSCYPO") * mm; //??
949    const double diameterHoleVoussFeet         = 2.0 * (*m_Abrt)[0]->getFloat("VOURCUTR") * mm;
950    const double thicknessTop    = (*m_Abrt)[0]->getFloat("VOUSPLOX") * mm; //???
951    const double thicknessCent   = (*m_Abrt)[0]->getFloat("VOUBZWTH") * mm;
952    const double thicknessBot    = (*m_Abrt)[0]->getFloat("VOUSPLIX") * mm; //???
953    const double thicknessRib         = (*m_Abrt)[0]->getFloat("VOURPYWI") * mm;
954    const double thicknessOutRibStd = (*m_Abrt)[0]->getFloat("VOUBLYWS") * mm;
955 //never used    const double ThetaAngleVoussFeet              = 15.0 * deg;
956 //!!
957 
958 //---------- BTVH -- BT Voussoir Head Parameters -----------------------------//
959 //                                                                            //
960                                                                               //
961    const double tHeightVHead         = (*m_Abrt)[0]->getFloat("CNBCOYEX") * mm;                    
962    const double tRMaxBTVoussHead     = (*m_Abrt)[0]->getFloat("CNBXMBRA") * mm;                    
963    const double xWidthMinVHead       = (*m_Abrt)[0]->getFloat("CNBCOXIN") * mm;                       
964    const double xWidthMaxVHead       = (*m_Abrt)[0]->getFloat("CNBCOXSU") * mm;                       
965    const double CnbEaXtp             = (*m_Abrt)[0]->getFloat("CNBEAXTP") * mm;
966 //never used    const double CnbEaYtp        = (*m_Abrt)[0]->getFloat("CNBEAYTP") * mm;
967    const double CnbEaXbt             = (*m_Abrt)[0]->getFloat("CNBEAXBT") * mm;
968 //never used    const double CnbEaYbt        = (*m_Abrt)[0]->getFloat("CNBEAYBT") * mm;
969    const double CnbEaCxi             = (*m_Abrt)[0]->getFloat("CNBEACXI") * mm;
970 //never used    const double CnbECYil        = (*m_Abrt)[0]->getFloat("CNBECYIL") * mm;
971    const double CnbECXol             = (*m_Abrt)[0]->getFloat("CNBECXOL") * mm;
972    const double CnbECXou             = (*m_Abrt)[0]->getFloat("CNBECXOU") * mm;
973 //never used    const double CnbECYol        = (*m_Abrt)[0]->getFloat("CNBECYOL") * mm;
974 //never used    const double CnbECYou        = (*m_Abrt)[0]->getFloat("CNBECYOU") * mm;
975 //never used    const double CnbECYiu        = (*m_Abrt)[0]->getFloat("CNBECYIU") * mm;
976    const double CnbEaCPl             = (*m_Abrt)[0]->getFloat("CNBEACPL") * mm;
977    const double CnbECZpo             = (*m_Abrt)[0]->getFloat("CNBECZPO") * mm;
978    const double CnbIECZp             = (*m_Abrt)[0]->getFloat("CNBIECZP") * mm;
979    const double CnbCaDma             = (*m_Abrt)[0]->getFloat("CNBCADMA") * mm;
980    const double CnbCaDmi             = (*m_Abrt)[0]->getFloat("CNBCADMI") * mm;
981    const double CnbCaZin             = (*m_Abrt)[0]->getFloat("CNBCAZIN") * mm;
982    const double CnbCaZex             = (*m_Abrt)[0]->getFloat("CNBCAZEX") * mm;
983    const double CnboxZex             = (*m_Abrt)[0]->getFloat("CNBOXZEX") * mm;
984 
985 //----------------------------------------------------------------------------
986   
987     
988   const double tolerance = 0.5 * mm;
989 
990   const double phi0 = M_PI/8;
991 
992   // array of phi values (centres of phi sectors)
993     const int nPhiSectors = 8;
994   double phiVal[nPhiSectors];
995   for ( int i = 0; i < nPhiSectors; i++ )  phiVal[i] = (2*i + 1) * phi0;  
996 
997   // auxiliary arrays to store phi values and positions for all volumes to be replicated 
998   const int nPosMax = maxDim;
999   double phiAux[ nPhiSectors * nPosMax ], posAux[ nPhiSectors * nPosMax ];
1000 
1001   const double XminEdge = (CnbEaXbt - xWidthMinVHead/2)/cos(phi0);
1002   const double xLengthWallEdge= (CnbEaXbt - CnbECXol)/cos(phi0);
1003   const double DeltaX = ( (tRMaxBTVoussHead - tHeightVHead) * sin(phi0) -
1004                        (xWidthMinVHead/2)*cos(phi0) - XminEdge - xLength/2 );
1005 //edge box
1006 //  GeoBox* sEdgeBox = new Box( xLengthWallEdge/2, yHeight/2, CnboxZex/2 );
1007   
1008 //  double distEdgeBox = xLength/2 + DeltaX + xLengthWallEdge/2)
1009   
1010 //    HepTransform3D trlEdgeBox[2];   
1011 //    trlEdgeBox[0] = HepTranslateX3D( -distEdgeBox );
1012 //    trlEdgeBox[1] = HepTranslateX3D( distEdgeBox );
1013                     
1014  //------------------------------------------------------------------------------------------
1015   // common Feet/Voussoir properties 
1016   const double yHeight = tWidthBTVouss;            // total y-height  
1017   const double yHeightCent   = yHeight - thicknessTop - thicknessBot;  // height of central plate
1018   const double radiusHole    = 0.5 * diameterHoleVoussFeet;  // radius of cutout hole
1019     
1020   // positioning variables
1021   // Feet Voussoir properties
1022 //  const double xLengthStd         = 2 * ( xPosVertex5StandFeet - thicknessInnerPlateStandFeet - xWidthConnVoussFeet);  // total length
1023   const double xLengthCentStd     = xLength - 2 * thicknessOutRibStd;  // length of central plate
1024 
1025   // positioning variables
1026   const double distOutRibStd    = 0.5 * ( xLength - 2 * thicknessOutRibStd - thicknessStepVouss);
1027   const double distOutRib  = 0.5 * ( xLength - thicknessOutRibStd - thicknessStepVouss );
1028   
1029   //------------------------------------------------------------------------------------------
1030   // central plate
1031   const GeoShape* sCentPlateStd = new GeoBox( xLengthCentStd/2 - thicknessStepVouss, yHeightCent/2, thicknessCent/2 );
1032 
1033   // top plate
1034   GeoBox* sTopPlateStd = new GeoBox( xLength/2 - tolerance, thicknessTop/2 - tolerance, zWidth/2 );
1035   
1036   // bottom plate
1037   GeoBox* sBotPlateStd = new GeoBox(xLength/2 - tolerance, thicknessBot/2 - tolerance, zWidth/2 );
1038 
1039 //edge box
1040   GeoBox* sEdgeBox = new GeoBox( xLengthWallEdge/2 - tolerance, yHeight/2, CnboxZex/2 );
1041   
1042   double distEdgeBox = xLength/2 + DeltaX + xLengthWallEdge/2;
1043       
1044   HepTransform3D trlEdgeBox[2];
1045   trlEdgeBox[0] = HepTranslate3D( -distEdgeBox, (thicknessTop - thicknessBot)/2, 0 );
1046   trlEdgeBox[1] = HepTranslate3D( distEdgeBox, (thicknessTop - thicknessBot)/2, 0 );
1047   
1048   //empty trapezoids  from top and bottom plates. 1-st trapezoid
1049   double xTopmin        = xLengthCentStd/2 - zWidthStepVouss/tan(alphaStep);
1050   double xTopmax        = xLengthCentStd/2;
1051   double yTopVouss      = yHeight;
1052   double zTopVouss      = zWidthStepVouss/2 + tolerance;
1053   GeoTrd* sEmptyTrap1   = new GeoTrd( xTopmin, xTopmax, yTopVouss, yTopVouss, zTopVouss );
1054   
1055   //empty trapezoid from top and bottom plates. 2-nd trapezoid
1056   double xBotmin        = xTopmax;
1057   double xBotmax        = xTopmin;
1058   double yBotVouss      = yHeight;
1059   double zBotVouss      = zTopVouss;
1060   GeoTrd* sEmptyTrap2   = new GeoTrd( xBotmin, xBotmax, yBotVouss, yBotVouss, zBotVouss );
1061   
1062   // inner rib pairs
1063   GeoBox* sInnerRib = new GeoBox( thicknessRib/2, yHeightCent/2 - thicknessStepVouss - tolerance, 
1064                                     zWidth/2 - zWidthStepVouss);  
1065   
1066   // outmost rib pair
1067   GeoBox* sOuterRibStd = new GeoBox( thicknessOutRibStd/2 + thicknessStepVouss/2 - tolerance, yHeightCent/2,
1068                                          zWidth/2 ); 
1069  
1070   //empty trapezoids at outer rib
1071   double xOutmin        = yHeightCent/2 - thicknessStepVouss;
1072   double xOutmax        = yHeightCent/2;
1073   double yOutVoussmax   = zWidth/2 - thicknessCent/2;
1074 //never used   double yOutVoussmin      = yOutVoussmax * (zWidthStepVouss - thicknessStepVouss * tan(alphaStep) )/zWidthStepVouss;
1075   double zOutVouss      = thicknessStepVouss/2 + tolerance;
1076   
1077   GeoTrd* sEmptyOuterRibTrap    = new GeoTrd( yOutVoussmax, yOutVoussmax, xOutmin, xOutmax, zOutVouss );
1078   
1079   // hole shape
1080   GeoTube* sHole = new GeoTube( 0, radiusHole, thicknessCent );
1081   
1082   //------------------------------------------------------------------------------------------
1083   // positions of shapes' centres relative to centre of central plate    
1084   HepTransform3D trlInnerRib[2];   // two inner rib pairs
1085   trlInnerRib[0] = HepTranslateX3D( -distInnerRib );
1086   trlInnerRib[1] = HepTranslateX3D( distInnerRib );
1087 
1088   
1089   HepTransform3D trlOutRibStd[2];   // outmost rib pair
1090   trlOutRibStd[0] = HepTranslateX3D( -distOutRib );
1091   trlOutRibStd[1] = HepTranslateX3D( distOutRib );
1092   
1093   HepTransform3D trlHole[2];   // hole pair
1094   trlHole[0] = HepTranslateX3D( -distHoleCentre );
1095   trlHole[1] = HepTranslateX3D( distHoleCentre );
1096    
1097   //position of empty shapes' centres relative to centre of central plate
1098   
1099   HepTransform3D trlEmptyOutRibTrap[4];
1100   trlEmptyOutRibTrap[0] = HepTranslate3D( distOutRibStd, 0, zWidth/2 + tolerance/2 ) * HepRotateY3D(3*M_PI/2);
1101   trlEmptyOutRibTrap[1] = HepTranslate3D( -distOutRibStd, 0, zWidth/2 + tolerance/2 ) * HepRotateY3D(M_PI/2);
1102   trlEmptyOutRibTrap[2] = HepTranslate3D( -distOutRibStd, 0, -zWidth/2 - tolerance/2 ) * HepRotateY3D(M_PI/2);
1103   trlEmptyOutRibTrap[3] = HepTranslate3D( distOutRibStd, 0, -zWidth/2 - tolerance/2 ) * HepRotateY3D(3*M_PI/2);
1104     
1105     // position of top/bottom plates and two empty trpezoids
1106     HepTransform3D trlTopPlate      = HepTranslate3D( 0, yHeightCent/2 + thicknessTop/2, 0 );
1107     HepTransform3D trlBotPlate      = HepTranslate3D( 0, -yHeightCent/2 - thicknessBot/2, 0 );
1108     
1109     HepTransform3D trlEmptyTrap1 = HepTranslate3D( 0, 0, zWidth/2 - zWidthStepVouss/2 );
1110     HepTransform3D trlEmptyTrap2 = HepTranslate3D( 0, 0, -zWidth/2 + zWidthStepVouss/2 );
1111 
1112   //------------------------------------------------------------------------------------------
1113   // add shapes to central plate, subtract holes
1114   
1115   // add inner ribs 
1116   for ( int i = 0; i < 2; i++ )
1117   {
1118      sCentPlateStd = &( sCentPlateStd->add( (*sInnerRib) << trlInnerRib[i] ) ); 
1119   }
1120   // add outer rib pair
1121   for ( int i = 0; i < 2; i++ )
1122   {
1123      sCentPlateStd = &( sCentPlateStd->add( (*sOuterRibStd) << trlOutRibStd[i] ) );
1124   } 
1125   // subtract holes
1126   for ( int i = 0; i < 2; i++ )
1127   {
1128      sCentPlateStd = &( sCentPlateStd->subtract( (*sHole) << trlHole[i] ) );
1129   } 
1130   //add EdgeWalls
1131   for ( int i = 0; i < 2; i++ )
1132     {
1133      sCentPlateStd = &( sCentPlateStd->add( (*sEdgeBox) << trlEdgeBox[i] ) );
1134     }
1135   //subtract shapes from outer rib
1136  
1137   for ( int i = 0; i < 4; i++ )
1138   {
1139      sCentPlateStd = &( sCentPlateStd->subtract( (*sEmptyOuterRibTrap) << trlEmptyOutRibTrap[i] ) );
1140   } 
1141   
1142   
1143   //-----------------------------------------------------------------------------------------------------------------------
1144   // Build Trapezoidal elements of inner rib
1145   //-----------------------------------------------------------------------------------------------------------------------
1146   
1147     //trapezoids at top and bottom edge of inner rib
1148 //     double xInnermin = thicknessRib/2;
1149 //     double xInnermax = thicknessRib/2 + thicknessStepVouss;
1150 //     double yInnerVouss       = zWidth/2 - zWidthStepVouss;
1151 //     double zInnerVouss       = thicknessStepVouss/2;
1152     
1153     //ss 12/02/2007 never user - memory never released const GeoShape* sTopInnerRibTrap = new GeoTrd( xInnermin, xInnermax, yInnerVouss, yInnerVouss, zInnerVouss );
1154     
1155     //ss 12/02/2007 never user - memory never released GeoTrd* sBotInnerRibTrap = new GeoTrd( xInnermax, xInnermin, yInnerVouss, yInnerVouss, zInnerVouss );
1156 
1157     HepTransform3D trlBot2TopInnerRibTrap = HepTranslate3D( 0, 0, -yHeightCent + thicknessStepVouss );
1158     
1159 //never used     const GeoShape& sTrapElementInnerRibVoussoir = sTopInnerRibTrap -> add( (*sBotInnerRibTrap) << trlBot2TopInnerRibTrap );
1160 
1161 
1162 //      GeoLogVol*  lTrapElementInnerRibVoussoir = new GeoLogVol( "TrapElementInnerRibVoussoir", 
1163 //                                              &sTrapElementInnerRibVoussoir, 
1164 //                                              theMaterialManager->getMaterial(voussMaterial) );
1165 //      GeoPhysVol* pTrapElementInnerRibVoussoir = new GeoPhysVol(lTrapElementInnerRibVoussoir);
1166 
1167 
1168   // centre of Voussoir central plate
1169   // (= reference point due to the assembly process chosen below)
1170   const double delta   = 0.5 * fabs( thicknessBot - thicknessTop );
1171   const double rCentre = 0.5 * ( yHeight + 2 * tRMinBTVouss ) - delta;
1172   const double centreXPos = rCentre * cos(phi0);
1173   const double centreYPos = -rCentre * sin(phi0);
1174 
1175   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1176   //!!! correction introduced for DC2, to void overlaps with BIR chambers (P03 layout) !!!
1177   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1178 //const double dc2Corr = zWidthSlimRib/2 - zWidth/2;
1179 
1180 //  GeoBox* sFatRib  = new GeoBox( thicknessFatRib/2, zWidth/2 + dc2Corr, tWidthCent/2 - tolerance );
1181 
1182   // positioning of Voussoir Head wrt. central plate   
1183 //HepTransform3D trlVoussHead = HepTranslate3D( (tRMaxBTVoussHead - tHeightVHead/2) * sin(phi0),
1184 //                    lengthHeadtoVous * cos(phi0),
1185 //                                              0 ); 
1186 
1187 const GeoShape& sBTVoussoir = sCentPlateStd->add( (*sTopPlateStd) << trlTopPlate ).
1188                                                add( (*sBotPlateStd) << trlBotPlate ).
1189                                                subtract( (*sEmptyTrap1) << trlEmptyTrap1 ).
1190                                                subtract( (*sEmptyTrap2) << trlEmptyTrap2 );
1191 
1192   GeoLogVol*  lBTVoussoir = new GeoLogVol( "BTVoussoir", 
1193                                              &sBTVoussoir, 
1194                                              theMaterialManager->getMaterial(voussoirMaterial) );
1195   GeoPhysVol* pBTVoussoir = new GeoPhysVol(lBTVoussoir);
1196 
1197 //edge of BTVoussoir
1198 
1199 //  const double XminEdge = (CnbEaXbt - xWidthMinVHead/2)/cos(phi0);
1200   const double XmaxEdge = (CnbEaXtp - xWidthMaxVHead/2)/cos(phi0);
1201   const double XminEmptyEdge1 = (CnbECXol - CnbEaCxi)/cos(phi0);
1202   const double XmaxEmptyEdge1 = (CnbECXou - CnbEaCxi)/cos(phi0);
1203   const double XminEmptyEdge2 = (CnbECXol - CnbEaCxi - CnbEaCPl)/cos(phi0);
1204   const double XmaxEmptyEdge2 = (CnbECXou - CnbEaCxi - CnbEaCPl)/cos(phi0);
1205                 
1206 //  const double DeltaX = ( (tRMaxBTVoussHead - tHeightVHead) * sin(phi0) - 
1207 //                      (xWidthMinVHead/2)*cos(phi0) - XminEdge - xLength/2 );
1208   
1209   
1210   const GeoShape* sEdgeVoussoir = new GeoTrd( XminEdge + DeltaX + xLength/2,
1211                                               XmaxEdge + DeltaX + xLength/2, 
1212                                               CnboxZex/2, CnboxZex/2,
1213                                               yHeight/2);
1214   
1215   GeoTrd* sEmptyTrap1EdgeBTVoussoir = new GeoTrd( XminEmptyEdge1 + DeltaX + xLength/2 + xLengthWallEdge,
1216                                                   XmaxEmptyEdge1 + DeltaX + xLength/2 + xLengthWallEdge,
1217                                                   CnboxZex/2, CnboxZex/2, yHeightCent/2 ); 
1218   GeoTrd* sEmptyTrap2EdgeBTVoussoir = new GeoTrd( XminEmptyEdge2 + DeltaX + xLength/2 + xLengthWallEdge,
1219                                                   XmaxEmptyEdge2 + DeltaX + xLength/2 + xLengthWallEdge,
1220                                                   CnboxZex/2, CnboxZex/2, yHeightCent/2 + tolerance );
1221     
1222   GeoBox* sEmptyBoxEdgeBTVoussoir = new GeoBox( DeltaX + xLength/2 + xLengthWallEdge, CnboxZex, yHeight);
1223   
1224   const GeoShape& sEdgeBTVoussoir = sEdgeVoussoir->subtract( (*sEmptyBoxEdgeBTVoussoir)
1225                                     << HepTranslate3D(0.,0.,0.) ).
1226                                     subtract( (*sEmptyTrap1EdgeBTVoussoir)
1227                                     << HepTranslate3D(0., CnbECZpo, -delta) ).
1228                                     subtract( (*sEmptyTrap1EdgeBTVoussoir)
1229                                     << HepTranslate3D(0., -CnbECZpo, -delta) ).
1230                                     subtract( (*sEmptyTrap2EdgeBTVoussoir)
1231                                     << HepTranslate3D(0., CnbIECZp, -delta) ).
1232                                     subtract( (*sEmptyTrap2EdgeBTVoussoir)
1233                                     << HepTranslate3D(0., -CnbIECZp, -delta) );
1234   
1235   GeoLogVol* lEdgeBTVoussoir = new GeoLogVol( "EdgeBTVoussoir", &sEdgeBTVoussoir,
1236                                              theMaterialManager->getMaterial(voussoirMaterial) );
1237   GeoPhysVol* pEdgeBTVoussoir = new GeoPhysVol(lEdgeBTVoussoir);
1238 
1239 
1240   // Voussoir Head properties 
1241 //never used   const double lengthHeadtoVous = (rCentre/cos(phi0)) + tHeightVHead/2 - tRMaxBTVoussHead;
1242 
1243   //------------------------------------------------------------------------------------------
1244   // build Voussoir Head as trapezoidal shape 
1245   const GeoShape* sVoussHead = new GeoTrd( xWidthMinVHead/2 - tolerance, xWidthMaxVHead/2 - tolerance,
1246                                    CnboxZex/2, CnboxZex/2, tHeightVHead/2 );  
1247   GeoCons* sEmptyHeadCons   = new GeoCons( 0., 0., CnbCaDma/2, CnbCaDmi/2, CnbCaZin/2 + tolerance,
1248                                          0., 2*M_PI );
1249   GeoTube* sEmptyHeadTube    = new GeoTube( 0., CnbCaDmi/2, CnbCaZex);
1250 
1251   const GeoShape& sHeadBTVoussoir = sVoussHead->subtract( (*sEmptyHeadCons)
1252                                  << HepTranslateZ3D(CnbCaZin/2 - tHeightVHead/2) ).
1253                                  subtract( (*sEmptyHeadTube) 
1254                                  << HepTranslateZ3D(-tHeightVHead/2) ).
1255                                  subtract( (*sEmptyTrap1EdgeBTVoussoir)
1256                                  << HepTranslate3D(centreXPos + delta*cos(phi0), CnbECZpo,
1257                                  -tRMaxBTVoussHead + tHeightVHead/2 + centreYPos - delta*sin(phi0))
1258                                  * HepRotateY3D(M_PI/2 + phi0) ).
1259                                  subtract( (*sEmptyTrap1EdgeBTVoussoir)
1260                                  << HepTranslate3D(centreXPos + delta*cos(phi0), -CnbECZpo,
1261                                  -tRMaxBTVoussHead + tHeightVHead/2 + centreYPos - delta*sin(phi0))
1262                                  * HepRotateY3D(M_PI/2 + phi0) ).
1263                                  subtract( (*sEmptyTrap1EdgeBTVoussoir)
1264                                  << HepTranslate3D(-centreXPos - delta*cos(phi0), CnbECZpo,
1265                                  -tRMaxBTVoussHead + tHeightVHead/2 + centreYPos - delta*sin(phi0))
1266                                  * HepRotateY3D(-M_PI/2 - phi0) ).
1267                                  subtract( (*sEmptyTrap1EdgeBTVoussoir)
1268                                  << HepTranslate3D(-centreXPos + delta*cos(phi0), -CnbECZpo,
1269                                  -tRMaxBTVoussHead + tHeightVHead/2 + centreYPos - delta*sin(phi0))
1270                                  * HepRotateY3D(-M_PI/2 - phi0) );
1271 
1272   GeoLogVol* lHeadBTVoussoir = new GeoLogVol( "HeadBTVoussoir", &sHeadBTVoussoir,
1273                                               theMaterialManager->getMaterial(voussoirMaterial) );
1274 
1275   GeoPhysVol* pHeadBTVoussoir = new GeoPhysVol(lHeadBTVoussoir);
1276 
1277   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1278   //!!! correction introduced for DC2, to avoid overlaps with BIR chambers (P03 layout) !!!
1279   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1280 //  GeoTrap* sBotSidePlate = new GeoTrap( tThicknessBot/2, thetaBotSide, 0,
1281 //                                        zWidth/2 + dc2Corr, dX1BotSide, dX1BotSide, 0,
1282 //                                        zWidth/2 + dc2Corr, dX3BotSide, dX3BotSide, 0 );
1283 
1284 
1285   
1286   // fill auxiliary arrays
1287   for ( int i = 0; i < nPhiSectors * nBTVouss; i++ )   
1288   {
1289      int ii = i % nPhiSectors, 
1290          jj = i / nPhiSectors;
1291      phiAux[i] = phiVal[ii]; 
1292      posAux[i] = zPosBTVouss[jj];
1293   }
1294 
1295   GENFUNCTION fVouss    = ArrayFunction( phiAux, phiAux + nPhiSectors * nBTVouss ); 
1296   GENFUNCTION fTrlVouss = ArrayFunction( posAux, posAux + nPhiSectors * nBTVouss ); 
1297   
1298   // central plate and Voussoir Head
1299   TRANSFUNCTION XFBTVoussoir = Pow( HepRotateZ3D(1.0), fVouss ) * 
1300                                  Pow( HepTranslateZ3D(1.0), fTrlVouss ) *
1301                                  HepTranslate3D( centreXPos, centreYPos, 0 ) *
1302                                  HepRotateZ3D(M_PI/2 - phi0);
1303                                  
1304   GeoSerialTransformer* sxBTVoussoir = new GeoSerialTransformer( pBTVoussoir,
1305                                                                    &XFBTVoussoir, 
1306                                                                    nPhiSectors * nBTVouss );
1307   container->add(sxBTVoussoir);
1308   
1309   TRANSFUNCTION XFEdgeBTVoussoir = Pow( HepRotateZ3D(1.0), fVouss ) *
1310                                     Pow( HepTranslateZ3D(1.0), fTrlVouss ) *
1311                                     HepTranslate3D( centreXPos + delta*cos(phi0),
1312                                     centreYPos - delta*sin(phi0), 0) * 
1313                                     HepRotateZ3D(M_PI/2 - phi0) *
1314                                     HepRotateX3D(M_PI/2);
1315   
1316   GeoSerialTransformer* sxEdgeBTVoussoir = new GeoSerialTransformer( pEdgeBTVoussoir,
1317                                                                      &XFEdgeBTVoussoir,
1318                                                                      nPhiSectors * nBTVouss );
1319   container->add(sxEdgeBTVoussoir);
1320 
1321   TRANSFUNCTION XFHeadBTVoussoir = Pow( HepRotateZ3D(1.0), fVouss ) *
1322                                     Pow( HepTranslateZ3D(1.0), fTrlVouss ) *
1323                                     HepTranslateX3D(tRMaxBTVoussHead - tHeightVHead/2) *
1324                                     HepRotateX3D(M_PI/2) * HepRotateY3D(M_PI/2);
1325 
1326   GeoSerialTransformer* sxHeadBTVoussoir = new GeoSerialTransformer( pHeadBTVoussoir,
1327                                                                      &XFHeadBTVoussoir,
1328                                                                      nPhiSectors * nBTVouss );
1329                                     
1330   container->add(sxHeadBTVoussoir);
1331 
1332 }
1333 
1334   void BarrelToroidBuilderRDB::buildBTStruts( GeoPhysVol* container ) 
1335 {
1336   ////////////////////////////////////////////////////////////////////////////////////////////
1337   //  Builds BT Struts (warm support structure)                                             //
1338   ////////////////////////////////////////////////////////////////////////////////////////////
1339 
1340   const StoredMaterialManager*  theMaterialManager;
1341   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
1342   {
1343     return;
1344   } 
1345 
1346   const int maxDim = 8; 
1347 
1348 
1349 //---------- BTVV -- BT Vacuum Vessel Parameters -----------------------------//
1350 //                                                                            //
1351                                                                               //
1352    const double rMinBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMIN") * mm;  //CryoRmin/101
1353    const double rMaxBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMAX") * mm;  //CryiRmax/102
1354    const double radiusBTVessel    = (*m_Abrt)[0]->getFloat("CRYORADI") * mm;  // 
1355 //                                                                            //
1356 //----------------------------------------------------------------------------//
1357 
1358 //---------- BTVS -- BT Voussoir Parameters ----------------------------------//
1359 //                                                                            //
1360                                                                               //
1361    int nBTVouss = 8;                                               //btvs->nvouss
1362    double zPosBTVouss[maxDim];                                                //
1363    for ( int i = 0; i < nBTVouss; i++ )  zPosBTVouss[i] = (*m_Abrt)[0]->getFloat("ZVOUSS",i) * mm; //
1364                                                                               //
1365 //----------------------------------------------------------------------------//
1366   
1367 //---------- BTST -- BT Strut Parameters -------------------------------------//
1368 //                                                                            //
1369                                                                               //
1370    const double tRMinBTStrut           = (*m_Abrt)[0]->getFloat("STRTRMIN") * mm;
1371    const double tRMaxBTStrut           = (*m_Abrt)[0]->getFloat("STRTRMAX") * mm;
1372    const double tThicknessBTStrut      = (*m_Abrt)[0]->getFloat("STRTRTHI") * mm;
1373    const double outerZWidthBTStrut     = (*m_Abrt)[0]->getFloat("STRTZWID") * mm;
1374    const double innerZThicknessBTStrut = (*m_Abrt)[0]->getFloat("STRTZTHI") * mm;
1375    const double lengthBTStrut          = (*m_Abrt)[0]->getFloat("STRTYLEN") * mm;
1376 //never used    const double dist2centre               = (*m_Abrt)[0]->getFloat("STRTRMOY") * mm;
1377    const double StrWRmax               = (*m_Abrt)[0]->getFloat("STRWRMAX") * mm;
1378    const double StrWYmax               = (*m_Abrt)[0]->getFloat("STRWYMAX") * mm;
1379    const double StrWRmin               = (*m_Abrt)[0]->getFloat("STRWRMIN") * mm;
1380    const double StrWYmin               = (*m_Abrt)[0]->getFloat("STRWYMIN") * mm;
1381    const double StrWZthi               = (*m_Abrt)[0]->getFloat("STRWZTHI") * mm;
1382    const double StrWYthi               = (*m_Abrt)[0]->getFloat("STRWYTHI") * mm;
1383    const double StrWZlen               = (*m_Abrt)[0]->getFloat("STRWZLEN") * mm;
1384    const double StrtYlnP               = (*m_Abrt)[0]->getFloat("STRTYLNP") * mm;
1385                                                                               //
1386    const std::string strutMaterial = getMaterial("Aluminium"); 
1387                                                                               //
1388 //----------------------------------------------------------------------------//
1389   
1390   const double tolerance = 0.5*mm;
1391 
1392   const double phi0 = M_PI/8;
1393 
1394   // array of phi values (centres of phi sectors)
1395   const int nPhiSectors = 8;
1396   double phiVal[nPhiSectors];
1397   for ( int i = 0; i < nPhiSectors; i++ )  phiVal[i] = ( 2*i + 1 ) * phi0;  
1398 
1399   // auxiliary arrays to store phi values and positions for all volumes to be replicated 
1400   const int nPosMax = maxDim;
1401   double phiAux[ nPhiSectors * nPosMax ], posAux[ nPhiSectors * nPosMax ];//, posAux1[ nPhiSectors * nPosMax ];
1402   
1403 
1404   // basic Strut properties 
1405   const double rOut             = radiusBTVessel; 
1406   const double halfLength       = lengthBTStrut/2;
1407   const double halfLengthbtubes = ( rMaxBTVessel - rOut ) * sin(phi0) - rOut;
1408   const double tWidth           = tRMaxBTStrut - tRMinBTStrut;
1409   const double tThickness       = tThicknessBTStrut;
1410   const double outerZWidth      = outerZWidthBTStrut;
1411   const double innerZThickness  = innerZThicknessBTStrut;
1412   
1413   // positioning variables for strut centre
1414   const double deltaR     = rMaxBTVessel - rMinBTVessel;
1415   const double rMean      = 0.5 * ( rMaxBTVessel + rMinBTVessel );
1416   const double centreXPos = deltaR/2 - rOut - ( halfLengthbtubes + rOut ) * sin(phi0); 
1417   const double centreYPos = -( halfLengthbtubes + rOut ) * cos(phi0);
1418 
1419 
1420   //------------------------------------------------------------------------------------------
1421   // Build Strut as bar with "I"-shape profile 
1422   // - let Struts touch the Vacuum Vessel as an approximation
1423   // - account for correct concave profile at sides touching the Vacuum Vessel 
1424   //------------------------------------------------------------------------------------------
1425 
1426 //  double alpha   = asin( tWidth/2/rOut );
1427 //never used   double sagitta = rOut - tWidth/2/tan(alpha);
1428 
1429   // inner vertical plate
1430   GeoBox* sVertBox = new GeoBox( halfLength - StrWYthi, tWidth/2 - tThickness - tolerance, innerZThickness/2 ); 
1431   
1432   // outer horizontal plates 
1433   GeoBox* sHorizBox = new GeoBox( halfLength - StrWYthi, tThickness/2, outerZWidth/2 ); 
1434 
1435   // empty tube to cut out the concave profiles at the sides
1436   //  GeoTube* sEmptyTube = new GeoTube( 0, rOut, outerZWidth/2 + tolerance );
1437   GeoBox* sEdgeBox = new GeoBox( StrWYthi/2, tWidth/2 + StrWZthi, StrWZlen/2 + StrWZthi );
1438   
1439 
1440   //------------------------------------------------------------------------------------------
1441   // unite horizontal and vertical plates, cut out left and right profiles
1442   const GeoShape& sStrut = sVertBox->add( (*sHorizBox) << HepTranslateY3D(  tWidth/2 - tThickness/2 ) ).
1443                                      add( (*sHorizBox) << HepTranslateY3D( -tWidth/2 + tThickness/2 ) ).
1444                                      add( (*sEdgeBox) << HepTranslateX3D( -halfLength + StrWYthi/2 ) ).
1445                                      add( (*sEdgeBox) << HepTranslateX3D( halfLength - StrWYthi/2 ) );
1446 //                                     subtract( (*sEmptyTube) << HepTranslateX3D( -halfLength - rOut ) ).
1447 //                                     subtract( (*sEmptyTube) << HepTranslateX3D(  halfLength + rOut ) );
1448 
1449   GeoLogVol*  lStrut = new GeoLogVol( "BTStrut", &sStrut, 
1450                                       theMaterialManager->getMaterial(strutMaterial) );
1451   GeoPhysVol* pStrut = new GeoPhysVol(lStrut);
1452 
1453   //----------------StrutWing-------------
1454 
1455   const double HeightTrap1 = StrWRmax - tRMaxBTStrut * cos(phi0) - sin(phi0) * StrtYlnP/2;
1456   const double HeightTrap2 = tRMaxBTStrut * cos(phi0) - tRMinBTStrut * cos(phi0);
1457   const double HeightTrap3 = tRMinBTStrut * cos(phi0) + sin(phi0) * StrtYlnP/2 - StrWRmin;
1458 
1459   GeoShape* sStrutWingTrap1  = new GeoTrd( tRMaxBTStrut * sin(phi0) - cos(phi0) * StrtYlnP/2,
1460                                            StrWYmax, StrWZthi/2 + StrWZlen/2, StrWZthi/2 + StrWZlen/2,
1461                                            HeightTrap1/2 );
1462   GeoTrd* sStrutWingTrap2    = new GeoTrd( tRMinBTStrut * sin(phi0) - cos(phi0) * StrtYlnP/2,
1463                                            tRMaxBTStrut * sin(phi0) - cos(phi0) * StrtYlnP/2,
1464                                            StrWZthi/2 + StrWZlen/2, StrWZthi/2 + StrWZlen/2,
1465                                            HeightTrap2/2 );
1466   GeoTrd* sStrutWingTrap3    = new GeoTrd( StrWYmin, tRMinBTStrut * sin(phi0) - 
1467                                            cos(phi0) * StrtYlnP/2,
1468                                            StrWZthi/2 + StrWZlen/2, StrWZthi/2 + StrWZlen/2, 
1469                                            HeightTrap3/2 );
1470   GeoTube* sEmptyStrutWingTube = new GeoTube( 0, rOut + tolerance, StrWZthi + StrWZlen+ tolerance );
1471   GeoBox* sEmptyStrutWingBox   = new GeoBox( tRMaxBTStrut * sin(phi0) - cos(phi0) * StrtYlnP/2 + 
1472                                              tolerance, StrWZlen/2 - StrWZthi/2, HeightTrap1 + 
1473                                              HeightTrap2 + HeightTrap3 + tolerance );
1474 
1475   const GeoShape& sWingStrut = sStrutWingTrap1->add( (*sStrutWingTrap2) << 
1476                                                 HepTranslateZ3D(-HeightTrap1/2 - HeightTrap2/2) ).
1477                                                 add( (*sStrutWingTrap3) << 
1478                                                 HepTranslateZ3D(-HeightTrap1/2 - HeightTrap2
1479                                                 - HeightTrap3/2) ).
1480                                                 subtract( (*sEmptyStrutWingTube) <<
1481                                                 HepTranslateZ3D( -HeightTrap1/2 )
1482                                                 * HepRotateX3D(M_PI/2) ).
1483                                                 subtract( (*sEmptyStrutWingBox) << 
1484                                                 HepTranslateZ3D( -(HeightTrap2 + HeightTrap3)/2 ) );
1485 
1486   GeoLogVol*  lWingStrut = new GeoLogVol( "BTWingStrut", &sWingStrut, 
1487                                       theMaterialManager->getMaterial(strutMaterial) );
1488   GeoPhysVol* pWingStrut = new GeoPhysVol(lWingStrut);
1489 
1490   //------------------------------------------------------------------------------------------
1491   // replicate and position
1492   for ( int i = 0; i < nPhiSectors * nBTVouss; i++ )  // fill auxiliary arrays  
1493   {
1494      int ii = i % nPhiSectors, 
1495          jj = i / nPhiSectors;
1496      phiAux[i] = phiVal[ii]; 
1497      posAux[i] = zPosBTVouss[jj];
1498    }
1499 
1500   GENFUNCTION fStrut    = ArrayFunction( phiAux, phiAux + nPhiSectors * nBTVouss ); 
1501   GENFUNCTION fTrlStrut = ArrayFunction( posAux, posAux + nPhiSectors * nBTVouss ); 
1502   GENFUNCTION fTrlWingStrut = ArrayFunction( posAux, posAux + nPhiSectors * nBTVouss );
1503   
1504   TRANSFUNCTION XFStrut = Pow( HepRotateZ3D(1.0), fStrut ) * 
1505                           Pow( HepTranslateZ3D(1.0), fTrlStrut ) *
1506                           HepTranslate3D( rMean + centreXPos, centreYPos, 0 ) *
1507                           HepRotateZ3D( -M_PI/2 - phi0 ); 
1508   TRANSFUNCTION XFWingStrut = Pow( HepRotateZ3D(1.0), fStrut ) *
1509                               Pow( HepTranslateZ3D(1.0), fTrlWingStrut ) *
1510                               HepTranslateX3D( rMean + deltaR/2 - radiusBTVessel + 
1511                                                HeightTrap1/2 ) *
1512                               HepRotateX3D(M_PI/2)* HepRotateY3D(M_PI/2);
1513 
1514 
1515   GeoSerialTransformer* sxStrut = new GeoSerialTransformer(pStrut, &XFStrut, nPhiSectors * nBTVouss);  GeoSerialTransformer* sxWingStrut = new GeoSerialTransformer(pWingStrut, &XFWingStrut,
1516                                                                nPhiSectors * nBTVouss);
1517 
1518   container->add(sxStrut);
1519   container->add(sxWingStrut);
1520   
1521 
1522 }
1523 
1524   void BarrelToroidBuilderRDB::buildBTCryoring( GeoPhysVol* container ) 
1525 {
1526   ////////////////////////////////////////////////////////////////////////////////////////////
1527   //  Builds BT Cryogenic Ring                                                              //
1528   ////////////////////////////////////////////////////////////////////////////////////////////
1529 
1530   const StoredMaterialManager*  theMaterialManager;
1531   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
1532   {
1533     return;
1534   } 
1535 
1536 //---------- BTVV -- BT Vacuum Vessel Parameters -----------------------------//
1537 //                                                                            //
1538                                                                               //
1539    const double rMinBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMIN") * mm;  //CryoRmin/101
1540    const double rMaxBTVessel      = (*m_Abrt)[0]->getFloat("CRYORMAX") * mm;  //CryiRmax/102
1541    const double radiusBTVessel    = (*m_Abrt)[0]->getFloat("CRYORADI") * mm;  // 
1542    const double thicknessBTVessel = (*m_Abrt)[0]->getFloat("CRYORADT") * mm;
1543 //                                                                            //
1544 //----------------------------------------------------------------------------//
1545 
1546 //---------- BTCR -- BT Cryogenic Ring Parameters ----------------------------//
1547 //                                                                            //
1548                                                                               //
1549    const double radiusBTCryoring = (*m_Abrt)[0]->getFloat("CRYRNGRA") * mm;   //CryoRingRadius/131 
1550  //never used   const double zPosBTCryoring   = (*m_Abrt)[0]->getFloat("CRYRNGZM") * mm;  //(XML bug)
1551    const double CryRngZm         = -1030.0 * mm;
1552                                                                    
1553 //est' esche i aluminium
1554    const std::string cryoringMaterial = getMaterial("Iron"); 
1555                                                                               //
1556 //----------------------------------------------------------------------------//
1557   
1558   
1559   // basic Cryoring(-tube) properties 
1560   const double phi0       = M_PI/8;
1561   const double rOut       = radiusBTVessel;
1562   const double rIn        = radiusBTCryoring - thicknessBTVessel; 
1563   const double halfLength = ( rMaxBTVessel - rOut ) * sin(phi0) - rOut;
1564   const double radius     = radiusBTCryoring;
1565 
1566   // positioning variables
1567   const double deltaR     = rMaxBTVessel - rMinBTVessel;
1568   const double rMean      = 0.5 * ( rMaxBTVessel + rMinBTVessel );
1569   const double centreXPos = deltaR/2 - rOut - ( halfLength + rOut ) * sin(phi0); 
1570   const double centreYPos = -( halfLength + rOut ) * cos(phi0);
1571   const double centreZPos = CryRngZm; 
1572 
1573 
1574   //------------------------------------------------------------------------------------------
1575   // Build Cryogenic Ring segment as tube 
1576   // - let Struts touch the Vacuum Vessel as an approximation
1577   // - account for correct concave profile at sides touching the Vacuum Vessel 
1578   //------------------------------------------------------------------------------------------
1579 
1580   double alpha   = asin( radius/rOut );
1581   double sagitta = rOut - radius/tan(alpha);
1582 
1583   GeoTube* sCryoringTube = new GeoTube( rIn, radius, halfLength + sagitta ); 
1584 
1585   // empty tube to cut out the concave profiles at the sides
1586   GeoTube* sEmptyTube = new GeoTube( 0, rOut, 2 * radius );
1587   
1588 
1589   //------------------------------------------------------------------------------------------
1590   // cut out left and right profiles
1591   
1592   // empty tubes have to be rotated
1593   HepTransform3D xfLeftTube  = HepTranslateZ3D( -halfLength - rOut ) * HepRotateX3D( M_PI/2 );
1594   HepTransform3D xfRightTube = HepTranslateZ3D(  halfLength + rOut ) * HepRotateX3D( M_PI/2 );
1595 
1596   const GeoShape& sCryoring = sCryoringTube->subtract( (*sEmptyTube) << xfLeftTube ).
1597                                              subtract( (*sEmptyTube) << xfRightTube );
1598 
1599   GeoLogVol*  lCryoring = new GeoLogVol( "BTCryoring", &sCryoring, theMaterialManager->getMaterial(cryoringMaterial) );
1600   GeoPhysVol* pCryoring = new GeoPhysVol(lCryoring);
1601 
1602 
1603   //------------------------------------------------------------------------------------------
1604   // replicate and position
1605   Variable iGen;
1606   GENFUNCTION f = ( 2*iGen + 1 ) * phi0;
1607   
1608   TRANSFUNCTION XFCryoring = Pow( HepRotateZ3D(1.0), f ) * 
1609                                   HepTranslate3D( rMean + centreXPos, centreYPos, centreZPos ) *
1610                                   HepRotateX3D( -M_PI/2 ) * HepRotateY3D(phi0); 
1611   GeoSerialTransformer* sxCryoring = new GeoSerialTransformer( pCryoring, &XFCryoring, 8 );
1612   container->add(sxCryoring);
1613 
1614 }
1615 
1616   void BarrelToroidBuilderRDB::buildRails( GeoPhysVol* container ) 
1617 {
1618   ////////////////////////////////////////////////////////////////////////////////////////////
1619   //  Builds Rails on Feet                                                                  //
1620   ////////////////////////////////////////////////////////////////////////////////////////////
1621 
1622   const StoredMaterialManager*  theMaterialManager;
1623   if ( StatusCode::SUCCESS != m_pDetStore->retrieve( theMaterialManager, "MATERIALS" ) ) 
1624   {
1625     return;
1626   } 
1627 
1628 //---------- FRLS -- Feet Rail Parameters ------------------------------------//
1629 //                                                                            //
1630                                                                               //
1631    const double zLengthCentralRail          = (*m_Rail)[0]->getFloat("CERZLENG") * mm;
1632    const double zLengthExtrRails            = (*m_Rail)[0]->getFloat("EXRZLENG") * mm;
1633    const double yMinRails                   = (*m_Rail)[0]->getFloat("YPOS")     * mm;
1634    const double yheightCentrPlateRail       = (*m_Rail)[0]->getFloat("CERTHIC2") * mm;
1635    const double xPosRails                   = (*m_Rail)[0]->getFloat("XPOS")     * mm;
1636    const double x2MinTopPlateRails          = (*m_Rail)[0]->getFloat("CERWID3I") * mm;
1637    const double x2MaxTopPlateRails          = (*m_Rail)[0]->getFloat("CERWID3O") * mm;
1638    const double xWidthBotPlateRails         = (*m_Rail)[0]->getFloat("CERWIDT1") * mm;
1639    const double thicknessHorizPlatesRails   = (*m_Rail)[0]->getFloat("CERTHIC1") * mm;
1640    const double thicknessVerticalPlateRails = (*m_Rail)[0]->getFloat("CERWIDT2") * mm;
1641    const double zLengthExtrRibsRail         = (*m_Rail)[0]->getFloat("CERRPEZL") * mm;
1642    const double zLengthCentralRibsRail      = (*m_Rail)[0]->getFloat("CERRPSZL") * mm;
1643    const double x2InPosEdgeRibsRails        = (*m_Rail)[0]->getFloat("CERRPIX1") * mm;
1644    const double x2OutPosEdgeRibsRails       = (*m_Rail)[0]->getFloat("CERRPOX1") * mm;
1645    const double angleInnerRibsCutOut        = 55.0                               * deg;
1646    const double angleOuterRibsCutOut        = 50.0                               * deg;
1647    const double zPosEdgeRibCentrRail1       = (*m_Rail)[0]->getFloat("CERRPE1Z") * mm;
1648    const double zPosEdgeRibCentrRail2       = (*m_Rail)[0]->getFloat("CERRPE2Z") * mm;
1649    const double zPosEdgeRibExtrRail1        = (*m_Rail)[0]->getFloat("EXRRPE1Z") * mm;
1650    const double zPosEdgeRibExtrRail2        = (*m_Rail)[0]->getFloat("EXRRPE2Z") * mm;
1651    
1652    const int nRibsCentralRail  = 11;
1653    const int nRibsExtrRails    = 17;
1654    
1655    double zPosRibsCentralRail[nRibsCentralRail];
1656    double zPosRibsExtrRails[nRibsExtrRails];
1657    
1658    for (int i = 0; i < nRibsCentralRail; i++ )
1659    {
1660    zPosRibsCentralRail[i] = (*m_Rail)[0]->getFloat("CERRPSZP",i) * mm;
1661    }
1662 
1663    for (int i = 0; i < nRibsExtrRails; i ++)
1664    {
1665    zPosRibsExtrRails[i]   = (*m_Rail)[0]->getFloat("EXRRPSZP",i) * mm;
1666    }                                                                             //
1667    const std::string railMaterial = getMaterial("Iron"); 
1668 //                                                                            //
1669 //----------------------------------------------------------------------------//
1670   
1671   const double tolerance = 0.5*mm;
1672 
1673   // general Ribs Rail properties
1674   const double xLengthRibRail   = x2InPosEdgeRibsRails + x2OutPosEdgeRibsRails
1675                                   + thicknessVerticalPlateRails;
1676   const double thicknessCentrRR = zLengthCentralRibsRail;
1677   const double thicknessExtrRR  = zLengthExtrRibsRail;
1678 
1679   // positioning ribs
1680   for (int i = 0; i <nRibsCentralRail; i++ )
1681   {
1682   zPosRibsCentralRail[i] = zPosRibsCentralRail[i] - zLengthCentralRibsRail/2;
1683   }
1684   for (int i = 0; i<nRibsExtrRails; i++ )
1685   {
1686   zPosRibsExtrRails[i] = zPosRibsExtrRails[i] - zLengthExtrRibsRail/2;
1687   }
1688   
1689   // general Rail properties
1690   const double zLength = zLengthCentralRail + 2 * zLengthExtrRails;
1691   const double yMin    = yMinRails; 
1692   const double yMax    = yMin + yheightCentrPlateRail + 2 * thicknessHorizPlatesRails; 
1693   
1694   // top+bottom plate properties
1695   const double thicknessPlates    = thicknessHorizPlatesRails;  
1696   
1697   // top plate properties 
1698   const double xMinTopPlate       = xPosRails - x2MinTopPlateRails; 
1699   const double xMaxTopPlate       = xPosRails + x2MaxTopPlateRails; 
1700   const double xWidthTopPlate     = xMaxTopPlate - xMinTopPlate; 
1701   
1702   // bottom plate properties 
1703   const double xWidthBotPlate     = xWidthBotPlateRails; 
1704   
1705   // vertical plate properties 
1706   const double thicknessVertPlate = thicknessVerticalPlateRails;  
1707   const double yHeightVertPlate   = yheightCentrPlateRail; 
1708   
1709   
1710   const double xylengthInRibBoxCut      = ( x2InPosEdgeRibsRails - xWidthBotPlate/2 + thicknessVertPlate/2) * cos(angleInnerRibsCutOut);
1711 
1712   const double xylengthOutRibBoxCut     = ( x2OutPosEdgeRibsRails - xWidthBotPlate/2 + thicknessVertPlate/2) * cos(angleOuterRibsCutOut);
1713 
1714   // positioning variables ...
1715   //  ... for vertical plate centre 
1716   const double xMean = xPosRails;   
1717   const double yMean = 0.5 * ( yMax + yMin );   
1718   
1719   //  ... for top plate centre (shifted in x)
1720   const double xShiftTopPlate = 0.5 * ( xMaxTopPlate + xMinTopPlate ) - xMean;   
1721 
1722 
1723   
1724   //------------------------------------------------------------------------------------------
1725   //  Build the Rails as union of the following box shapes:
1726   //   - central, vertical plate
1727   //   - horizontal bottom plate, aligned in x with central plate
1728   //   - horizontal top plate, shifted outside in x wrt. the central plate
1729   //------------------------------------------------------------------------------------------
1730 
1731   // central, vertical plate
1732   const GeoShape* sCentPlate = new GeoBox( thicknessVertPlate/2, yHeightVertPlate/2 - tolerance, zLength/2 ); 
1733   
1734   // bottom plate
1735   GeoBox* sBotPlate  = new GeoBox( xWidthBotPlate/2, thicknessPlates/2 - tolerance, zLength/2 ); 
1736   
1737   // top plate
1738   GeoBox* sTopPlate  = new GeoBox( xWidthTopPlate/2, thicknessPlates/2 - tolerance, zLength/2 ); 
1739   
1740   // central ribs
1741   const GeoShape* sCentrRib  = new GeoBox( xLengthRibRail/2, yHeightVertPlate/2 - tolerance, thicknessCentrRR/2 );
1742   
1743   // extr ribs
1744   const GeoShape* sExtrRib   = new GeoBox(xLengthRibRail/2, yHeightVertPlate/2 - tolerance, thicknessExtrRR/2 );
1745   
1746   //empty central box from ribs
1747   GeoBox* sEmptyCBoxRib   = new GeoBox(thicknessVertPlate/2 + tolerance,
1748                                        yHeightVertPlate/2, thicknessCentrRR );
1749   
1750   HepTransform3D trlEmptyCBoxRib = HepTranslateX3D( (x2OutPosEdgeRibsRails - x2InPosEdgeRibsRails)/2 );
1751 
1752   //empty inner&outer boxes from ribs
1753   GeoBox* sEmptyInBoxRib  = new GeoBox(xylengthInRibBoxCut, 2*xylengthInRibBoxCut, thicknessCentrRR );
1754   GeoBox* sEmptyOutBoxRib = new GeoBox(xylengthOutRibBoxCut, 2*xylengthOutRibBoxCut, thicknessCentrRR );
1755     
1756   HepTransform3D trlEmptyInBoxRib = HepTranslate3D( xLengthRibRail/2, -yHeightVertPlate/2, 0. )
1757                                                     * HepRotateZ3D(-angleInnerRibsCutOut);
1758 
1759   HepTransform3D trlEmptyOutBoxRib = HepTranslate3D( -xLengthRibRail/2, -yHeightVertPlate/2, 0. )
1760                                                      * HepRotateZ3D(angleOuterRibsCutOut);
1761   // build rails ribs
1762   sCentrRib = &( sCentrRib->subtract( (*sEmptyCBoxRib) << trlEmptyCBoxRib ).
1763                             subtract( (*sEmptyInBoxRib) << trlEmptyInBoxRib ).
1764                             subtract( (*sEmptyOutBoxRib) << trlEmptyOutBoxRib ) );
1765 
1766   const double deltax = -(x2OutPosEdgeRibsRails - x2InPosEdgeRibsRails)/2;
1767 
1768   sExtrRib  = &( sExtrRib->subtract( (*sEmptyCBoxRib) << trlEmptyCBoxRib ).
1769                            subtract( (*sEmptyInBoxRib) << trlEmptyInBoxRib ).
1770                            subtract( (*sEmptyOutBoxRib) << trlEmptyOutBoxRib ) );
1771 
1772   for (int i = 0; i < nRibsCentralRail; i++ )
1773   {
1774   sCentPlate = & (sCentPlate->add( (*sCentrRib) << HepTranslate3D( deltax, 0., zPosRibsCentralRail[i]
1775                                                                   - zLengthCentralRail/2) ) );
1776   }
1777   
1778   for (int i = 0; i < nRibsExtrRails; i++ )
1779   {
1780   sCentPlate = &( sCentPlate->add( (*sExtrRib) << HepTranslate3D( deltax, 0., zPosRibsExtrRails[i]
1781                                                                 - zLengthCentralRail/2
1782                                                                 - zLengthExtrRails) ) );
1783   sCentPlate = &( sCentPlate->add( (*sExtrRib) << HepTranslate3D( deltax, 0., zPosRibsExtrRails[i]
1784                                                                  + zLengthCentralRail/2) ) );
1785   }
1786 // add ribs to the edge of central rail part and 2 extremity rail parts
1787 
1788   sCentPlate = &( sCentPlate->add( (*sExtrRib) << HepTranslate3D( deltax, 0., zPosEdgeRibCentrRail1 ) ).
1789                               add( (*sExtrRib) << HepTranslate3D( deltax, 0., zPosEdgeRibCentrRail2 ) ).
1790                               add( (*sExtrRib) << HepTranslate3D( deltax, 0., -zLengthCentralRail/2 - zLengthExtrRails/2 + zPosEdgeRibExtrRail1 ) ).
1791                               add( (*sExtrRib) << HepTranslate3D( deltax, 0., -zLengthCentralRail/2 - zLengthExtrRails/2 + zPosEdgeRibExtrRail2 ) ).
1792                               add( (*sExtrRib) << HepTranslate3D( deltax, 0., zLengthCentralRail/2 + zLengthExtrRails/2 + zPosEdgeRibExtrRail1 ) ).
1793                               add( (*sExtrRib) << HepTranslate3D( deltax, 0., zLengthCentralRail/2 + zLengthExtrRails/2 + zPosEdgeRibExtrRail2 ) ) );
1794                                                      
1795   //------------------------------------------------------------------------------------------
1796   // add bottom/top plates to central plate
1797   double yPosPlates = 0.5 * ( yHeightVertPlate + thicknessPlates );
1798   const GeoShape& sRail = sCentPlate->add( (*sBotPlate) << HepTranslateY3D( -yPosPlates ) ).
1799                                       add( (*sTopPlate) << HepTranslate3D( -xShiftTopPlate, yPosPlates, 0 ) );
1800 
1801   GeoLogVol*  lRail = new GeoLogVol( "Rail", &sRail, theMaterialManager->getMaterial(railMaterial) );
1802   GeoPhysVol* pRail = new GeoPhysVol(lRail);
1803 
1804 
1805   //------------------------------------------------------------------------------------------
1806   // replicate and position
1807   //------------------------------------------------------------------------------------------
1808   // array of x-positions
1809   const int nRails = 2;
1810   double xPos[nRails] = { -xMean, xMean };
1811 
1812   // array of rotation angles (flipping around y-axis)
1813   double rotAngle[nRails] = { 0, M_PI };
1814 
1815   GENFUNCTION fRot = ArrayFunction( rotAngle, rotAngle + nRails );
1816   GENFUNCTION fTrl = ArrayFunction( xPos,     xPos     + nRails );
1817 
1818   TRANSFUNCTION XFRail = Pow( HepTranslateX3D(1.0), fTrl ) * 
1819                          HepTranslateY3D(yMean) *             // constant displacement     
1820                          Pow( HepRotateY3D(1.0), fRot );
1821   
1822   GeoSerialTransformer* sxRail = new GeoSerialTransformer( pRail, &XFRail, nRails );
1823   container->add(sxRail);
1824 
1825 } 
1826 
1827 
1828 
1829 std::string BarrelToroidBuilderRDB::getMaterial( std::string materialName ) 
1830 {
1831   MsgStream log(Athena::getMessageSvc(), "MuonGeoModel");
1832   if ( materialName == "Alum" )
1833   {
1834     return "std::Aluminium";
1835   }
1836   else if ( materialName == "Iron" )  
1837   {
1838     return "std::Iron";
1839   }  
1840   else if ( materialName == "Fe50"  ||  materialName == "Al67" )   
1841   {
1842     return "toro::" + materialName;
1843   }
1844   else 
1845   {
1846     log  <<  "Barrel ToroidBuilderRDB::getMaterial: material "  <<  materialName  <<  " not defined! "  
1847          <<  " Take Aluminium instead." 
1848          <<  endreq;
1849     return "std::Aluminium";           
1850   }  
1851 }  
1852 
1853 
1854 } // namespace MuonGM
1855
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine. Valid HTML 4.01!