/DetectorDescription/Parser/src/DDDividedTrd.cc

  1//
  2// ********************************************************************
  3// 25.04.04 - M. Case ddd-ize G4ParameterisationTrd*
  4// ********************************************************************
  5
  6#include "DetectorDescription/Parser/src/DDDividedTrd.h"
  7#include "DetectorDescription/Parser/src/DDXMLElement.h"
  8
  9#include "DetectorDescription/Core/interface/DDLogicalPart.h"
 10#include "DetectorDescription/Core/interface/DDName.h"
 11#include "DetectorDescription/Core/interface/DDAxes.h"
 12#include "DetectorDescription/Core/interface/DDSolid.h"
 13#include "DetectorDescription/Core/interface/DDMaterial.h"
 14
 15#include "DetectorDescription/Base/interface/DDdebug.h"
 16
 17#include "CLHEP/Units/GlobalSystemOfUnits.h"
 18
 19#include <cmath>
 20#include <cstdlib>
 21
 22DDDividedTrdX::DDDividedTrdX( const DDDivision& div, DDCompactView* cpv )
 23  : DDDividedGeometryObject(div,cpv)
 24{
 25  checkParametersValidity();
 26  setType( "DivisionTrdX" );
 27  DDTrap mtrd = (DDTrap)( div_.parent().solid() );
 28
 29  if ( divisionType_ == DivWIDTH )
 30  {
 31    compNDiv_ = calculateNDiv( 2 * mtrd.x1(), div_.width(), div_.offset() );
 32  }
 33  else if( divisionType_ == DivNDIV )
 34  {
 35    compWidth_ = calculateWidth( 2*mtrd.x1(), div_.nReplicas(), div_.offset() );
 36  }
 37
 38  DCOUT_V ('P', " DDDividedTrdX - ## divisions " << compNDiv_ << " = " << div_.nReplicas() << "\n Offset " << div_.offset() << "\n Width " << compWidth_ << " = " << div_.width());
 39}
 40
 41DDDividedTrdX::~DDDividedTrdX( void )
 42{}
 43
 44double
 45DDDividedTrdX::getMaxParameter( void ) const
 46{
 47  DDTrap mtrd = (DDTrap)(div_.parent().solid());
 48  return 2 * mtrd.x1();
 49}
 50
 51DDTranslation
 52DDDividedTrdX::makeDDTranslation( const int copyNo ) const
 53{
 54  DDTrap mtrd = (DDTrap)(div_.parent().solid());
 55  double mdx = mtrd.x1();
 56
 57
 58  //----- translation 
 59  double posi = -mdx + div_.offset() + (copyNo+0.5)*compWidth_;
 60
 61  DCOUT_V ('P', " DDDividedTrdX: " << copyNo << "\n Position: x=" << posi << "  Axis= " << DDAxesNames::name(div_.axis()) << "\n");
 62
 63  if( div_.axis() == x )
 64  {
 65    return DDTranslation(posi, 0.0, 0.0);
 66  }
 67  else
 68  { 
 69    std::string s = "ERROR - DDDividedTrdX::makeDDTranslation()";
 70    s += "\n        Axis is along ";
 71    s += DDAxesNames::name(div_.axis());
 72    s += " !\n" ;
 73    s += "DDDividedTrdX::makeDDTranslation()";
 74    s += " IllegalConstruct: Only axes along x are allowed !";
 75    throw cms::Exception("DDException") << s;
 76  }
 77  
 78  return DDTranslation();
 79}
 80
 81DDRotation
 82DDDividedTrdX::makeDDRotation( const int copyNo ) const
 83{
 84  return DDRotation();
 85}
 86
 87DDLogicalPart
 88DDDividedTrdX::makeDDLogicalPart( const int copyNo ) const
 89{
 90  DDTrap mtrd = (DDTrap)(div_.parent().solid());
 91  DDMaterial usemat = div_.parent().material();
 92
 93  double pDy1 = mtrd.y1(); //GetYHalfLength1();
 94  double pDy2 = mtrd.y2(); //->GetYHalfLength2();
 95  double pDz = mtrd.halfZ(); //->GetZHalfLength();
 96  double pDx = compWidth_/2.;
 97 
 98  //trd.SetAllParameters ( pDx, pDx, pDy1, pDy2, pDz );
 99
100  DDName solname(div_.parent().ddname().name() + "_DIVCHILD" 
101		 , div_.parent().ddname().ns());
102  DDSolid dsol(solname);
103  DDLogicalPart ddlp(solname);
104  if (!dsol.isDefined().second)
105  {
106    dsol = DDSolidFactory::trap(solname
107				, pDz
108				, 0.*deg
109				, 0.*deg
110				, pDy1
111				, pDx
112				, pDx
113				, 0.*deg
114				, pDy2
115				, pDx
116				, pDx
117				, 0.*deg);
118    ddlp = DDLogicalPart(solname, usemat, dsol);
119  }
120  DCOUT_V ('P', "DDDividedTrdX::makeDDLogicalPart lp = " << ddlp);
121  return ddlp;
122}
123
124void
125DDDividedTrdX::checkParametersValidity( void )
126{
127  DDDividedGeometryObject::checkParametersValidity();
128  
129  DDTrap mtrd = (DDTrap)(div_.parent().solid());
130
131  double mpDx1 = mtrd.x1(); //->GetXHalfLength1();
132  double mpDx2 = mtrd.x2(); //->GetXHalfLength2();
133  double mpDx3 = mtrd.x3(); 
134  double mpDx4 = mtrd.x4();
135  double mpTheta = mtrd.theta();
136  double mpPhi = mtrd.phi();
137  double mpAlpha1 = mtrd.alpha1();  
138  double mpAlpha2 = mtrd.alpha2();
139  //    double mpDy1 = mtrd.y1();
140  //    double mpDy2 = mtrd.y2();
141  //   double x1Tol = mpDx1 - mpDx2;
142  //   if (x1Tol < 0.0) x1Tol = x1Tol * -1.0;
143
144  if ( fabs(mpDx1 - mpDx2) > tolerance()  || fabs(mpDx3 - mpDx4) > tolerance()
145       || fabs(mpDx1 - mpDx4) > tolerance())
146  {
147    std::string s = "ERROR - DDDividedTrdX::checkParametersValidity()";
148    s+= "\n        Making a division of a TRD along axis X,";
149    s+= "\n        while the X half lengths are not equal,";
150    s+= "\n        is not (yet) supported. It will result";
151    s+= "\n        in non-equal division solids.";
152    throw cms::Exception("DDException") << s;
153  }
154  //    if (fabs(mpDy1 - mpDy2) > tolerance())
155  //      {
156  //        std::string s = "ERROR - DDDividedTrdX::checkParametersValidity()";
157  //        s+= "\n        Making a division of a TRD along axis X,";
158  //        s+= "\n        while the Y half lengths are not equal,";
159  //        s+= "\n        is not (yet) supported. It will result";
160  //        s+= "\n        in non-equal division solids.";
161  //        throw cms::Exception("DDException") << s;
162  //      }
163  // mec:  we only have traps, not trds in DDD, so I added this check
164  // to make sure it is only a trd (I think! :-))
165  if (mpAlpha1 != 0.*deg || mpAlpha2 != 0.*deg || mpTheta != 0.*deg || mpPhi != 0.*deg)
166  {
167    std::string s = "ERROR - DDDividedTrdX::checkParametersValidity()";
168    s+= "\n        Making a division of a TRD along axis X,";
169    s+= "\n        while the theta, phi and aplhpa2 are not zero,";
170    s+= "\n        is not (yet) supported. It will result";
171    s+= "\n        in non-equal division solids.";
172    throw cms::Exception("DDException") << s;
173  }
174}
175
176DDDividedTrdY::DDDividedTrdY( const DDDivision& div, DDCompactView* cpv )
177  : DDDividedGeometryObject( div, cpv )
178{
179  checkParametersValidity();
180  setType( "DivisionTrdY" );
181  DDTrap mtrd = (DDTrap)(div_.parent().solid());
182
183  if( divisionType_ == DivWIDTH )
184  {
185    compNDiv_ = calculateNDiv( 2 * mtrd.y1(), div_.width(), div_.offset() );
186  }
187  else if( divisionType_ == DivNDIV )
188  {
189    compWidth_ = calculateWidth( 2 * mtrd.y1(), div_.nReplicas(), div_.offset() );
190  }
191
192  DCOUT_V ('P', " DDDividedTrdY no divisions " << compNDiv_ << " = " << div_.nReplicas() << "\n Offset " << div_.offset() << "\n width " << compWidth_ << " = " << div_.width() << std::endl);  
193}
194
195DDDividedTrdY::~DDDividedTrdY( void )
196{}
197
198double
199DDDividedTrdY::getMaxParameter( void ) const
200{
201  DDTrap mtrd = (DDTrap)(div_.parent().solid());
202  return 2 * mtrd.y1(); 
203}
204
205DDTranslation
206DDDividedTrdY::makeDDTranslation( const int copyNo ) const
207{
208  DDTrap mtrd = (DDTrap)(div_.parent().solid() );
209  double mdy = mtrd.y1();
210
211  //----- translation 
212  double posi = -mdy + div_.offset() + (copyNo+0.5)*compWidth_;
213
214  DCOUT_V ('P', " DDDividedTrdY: " << copyNo << "\n Position: y=" << posi << "  Axis= " << DDAxesNames::name(div_.axis()) << "\n");
215
216  if( div_.axis() == y )
217  {
218    return DDTranslation(0.0, posi, 0.0);
219  }
220  else
221  { 
222    std::string s = "ERROR - DDDividedTrdY::makeDDTranslation()";
223    s += "\n        Axis is along ";
224    s += DDAxesNames::name(div_.axis());
225    s += " !\n" ;
226    s += "DDDividedTrdY::makeDDTranslation()";
227    s += " IllegalConstruct: Only axes along y are allowed !";
228    throw cms::Exception("DDException") << s;
229  }
230  return DDTranslation();
231}
232
233DDRotation
234DDDividedTrdY::makeDDRotation( const int copyNo ) const
235{
236  return DDRotation();
237}
238
239DDLogicalPart
240DDDividedTrdY::makeDDLogicalPart( const int copyNo ) const
241{
242  //---- The division along Y of a Trd will result a Trd, only 
243  //--- if Y at -Z and +Z are equal, else use the G4Trap version
244  DDTrap mtrd = (DDTrap)(div_.parent().solid());
245  DDMaterial usemat = div_.parent().material();
246  
247  double pDx1 = mtrd.x1(); //->GetXHalfLength1() at Y+;
248  double pDx2 = mtrd.x2(); //->GetXHalfLength2() at Y+;
249  double pDx3 = mtrd.x3(); //->GetXHalfLength1() at Y-;
250  double pDx4 = mtrd.x4(); //->GetXHalfLength2() at Y-;
251  double pDz = mtrd.halfZ(); //->GetZHalfLength();
252  double pDy = compWidth_/2.;
253 
254  //trd.SetAllParameters ( pDx1, pDx2, pDy, pDy, pDz );
255  DDName solname(div_.name() );
256  DDSolid  dsol(solname);
257  DDLogicalPart ddlp(solname);
258  if (!dsol.isDefined().second)
259  {
260    dsol = DDSolidFactory::trap(solname
261				, pDz
262				, 0.*deg
263				, 0.*deg
264				, pDy
265				, pDx1
266				, pDx2
267				, 0*deg
268				, pDy
269				, pDx3
270				, pDx4
271				, 0.*deg);
272    DDLogicalPart ddlp(solname,  usemat, dsol);
273  }
274  DCOUT_V ('P', "DDDividedTrdY::makeDDLogicalPart lp = " << ddlp);
275  return ddlp;
276}
277
278void
279DDDividedTrdY::checkParametersValidity( void )
280{
281  DDDividedGeometryObject::checkParametersValidity();
282
283  DDTrap mtrd = (DDTrap)(div_.parent().solid());
284
285  double mpDy1 = mtrd.y1(); //->GetYHalfLength1();
286  double mpDy2 = mtrd.y2(); //->GetYHalfLength2();
287  double mpTheta = mtrd.theta();
288  double mpPhi = mtrd.phi();
289  double mpAlpha1 = mtrd.alpha1();
290  double mpAlpha2 = mtrd.alpha2();
291
292  if( fabs(mpDy1 - mpDy2) > tolerance() )
293  {
294    std::string s= "ERROR - DDDividedTrdY::checkParametersValidity()";
295    s += "\n        Making a division of a TRD along axis Y while";
296    s += "\n        the Y half lengths are not equal is not (yet)";
297    s += "\n        supported. It will result in non-equal";
298    s += "\n        division solids.";
299    throw cms::Exception("DDException") << s;
300  }
301  // mec:  we only have traps, not trds in DDD, so I added this check
302  // to make sure it is only a trd (I think! :-))
303  if (mpAlpha1 != 0.*deg || mpAlpha2 != 0.*deg || mpTheta != 0.*deg || mpPhi != 0.*deg)
304  {
305    std::string s = "ERROR - DDDividedTrdY::checkParametersValidity()";
306    s+= "\n        Making a division of a TRD along axis X,";
307    s+= "\n        while the theta, phi and aplhpa2 are not zero,";
308    s+= "\n        is not (yet) supported. It will result";
309    s+= "\n        in non-equal division solids.";
310    throw cms::Exception("DDException") << s;
311  }
312}
313
314DDDividedTrdZ::DDDividedTrdZ( const DDDivision& div, DDCompactView* cpv )
315  : DDDividedGeometryObject( div, cpv )
316{ 
317  checkParametersValidity();
318  setType( "DivTrdZ" );
319  DDTrap mtrd = (DDTrap)(div_.parent().solid());
320
321  if ( divisionType_ == DivWIDTH )
322  {
323    compNDiv_ = calculateNDiv( 2*mtrd.halfZ(), div_.width(), div_.offset() );
324  }
325  else if( divisionType_ == DivNDIV )
326  {
327    compWidth_ = calculateWidth( 2*mtrd.halfZ(), div_.nReplicas(), div_.offset() );
328  }
329  DCOUT_V ('P', " DDDividedTrdY no divisions " << compNDiv_ << " = " << div_.nReplicas() << "\n Offset " << div_.offset() << "\n width " << compWidth_ << " = " << div_.width() << std::endl);
330}
331
332DDDividedTrdZ::~DDDividedTrdZ( void )
333{}
334
335double
336DDDividedTrdZ::getMaxParameter( void ) const
337{
338  DDTrap mtrd = (DDTrap)(div_.parent().solid());
339  return 2 * mtrd.halfZ();
340}
341
342DDTranslation
343DDDividedTrdZ::makeDDTranslation( const int copyNo ) const
344{
345  DDTrap mtrd = (DDTrap)(div_.parent().solid() );
346  double mdz = mtrd.halfZ();
347
348  //----- translation 
349  double posi = -mdz + div_.offset() + (copyNo+0.5)*compWidth_;
350
351  DCOUT_V ('P', " DDDividedTrdZ: " << copyNo << "\n Position: z=" << posi << "  Axis= " << DDAxesNames::name(div_.axis()) << "\n");
352
353  if( div_.axis() == z )
354  {
355    return DDTranslation(0.0, 0.0, posi);
356  }
357  else
358  { 
359    std::string s = "ERROR - DDDividedTrdZ::makeDDTranslation()";
360    s += "\n        Axis is along ";
361    s += DDAxesNames::name(div_.axis());
362    s += " !\n" ;
363    s += "DDDividedTrdY::makeDDTranslation()";
364    s += " IllegalConstruct: Only axes along z are allowed !";
365    throw cms::Exception("DDException") << s;
366
367  }
368  return DDTranslation();
369}
370
371DDRotation
372DDDividedTrdZ::makeDDRotation( const int copyNo ) const
373{
374  return DDRotation();
375}
376
377DDLogicalPart
378DDDividedTrdZ::makeDDLogicalPart ( const int copyNo ) const
379{
380  //---- The division along Z of a Trd will result a Trd
381  DDTrap mtrd = (DDTrap)(div_.parent().solid());
382  DDMaterial usemat = div_.parent().material();
383
384  double pDx1 = mtrd.x1(); //->GetXHalfLength1();
385  //(mtrd->GetXHalfLength2() - mtrd->GetXHalfLength1() );
386  double DDx = (mtrd.x2() - mtrd.x1() );
387  double pDy1 = mtrd.y1(); // ->GetYHalfLength1();
388  //(mtrd->GetYHalfLength2() - mtrd->GetYHalfLength1() );
389  double DDy = (mtrd.y2() - mtrd.y1() );
390  double pDz = compWidth_/2.;
391  double zLength = 2*mtrd.halfZ(); //->GetZHalfLength();
392 
393  //    trd.SetAllParameters ( pDx1+DDx*(div_.offset()+copyNo*compWidth_)/zLength,
394  //                           pDx1+DDx*(div_.offset()+(copyNo+1)*compWidth_)/zLength, 
395  //                           pDy1+DDy*(div_.offset()+copyNo*compWidth_)/zLength,
396  //                           pDy1+DDy*(div_.offset()+(copyNo+1)*compWidth_)/zLength, pDz );
397
398  DDName solname(div_.parent().ddname().name() + "_DIVCHILD" 
399		 + DDXMLElement::itostr(copyNo)
400		 , div_.parent().ddname().ns());
401  DDSolid  dsol = 
402    DDSolidFactory::trap(solname
403			 , pDz
404			 , 0.*deg
405			 , 0.*deg
406			 , pDy1+DDy*(div_.offset()+copyNo*compWidth_)/zLength
407			 , pDx1+DDx*(div_.offset()+copyNo*compWidth_)/zLength
408			 , pDx1+DDx*(div_.offset()+copyNo*compWidth_)/zLength
409			 , 0.*deg
410			 , pDy1+DDy*(div_.offset()+(copyNo+1)*compWidth_)/zLength
411			 , pDx1+DDx*(div_.offset()+(copyNo+1)*compWidth_)/zLength
412			 , pDx1+DDx*(div_.offset()+(copyNo+1)*compWidth_)/zLength
413			 , 0*deg
414      );
415
416  DDLogicalPart ddlp(solname, usemat, dsol);
417  DCOUT_V ('P', "DDDividedTrdZ::makeDDLogicalPart lp = " << ddlp);
418  return ddlp;
419}
420
421void
422DDDividedTrdZ::checkParametersValidity( void )
423{
424  DDDividedGeometryObject::checkParametersValidity();
425
426  DDTrap mtrd = (DDTrap)(div_.parent().solid());
427
428  double mpTheta = mtrd.theta();
429  double mpPhi = mtrd.phi();
430  double mpAlpha1 = mtrd.alpha1();
431  double mpAlpha2 = mtrd.alpha2();
432
433  // mec:  we only have traps, not trds in DDD, so I added this check
434  // to make sure it is only a trd (I think! :-))
435  if (mpAlpha1 != 0.*deg || mpAlpha2 != 0.*deg || mpTheta != 0.*deg || mpPhi != 0.*deg)
436  {
437    std::string s = "ERROR - DDDividedTrdZ::checkParametersValidity()";
438    s+= "\n        Making a division of a TRD along axis X,";
439    s+= "\n        while the theta, phi and aplhpa2 are not zero,";
440    s+= "\n        is not (yet) supported. It will result";
441    s+= "\n        in non-equal division solids.";
442    throw cms::Exception("DDException") << s;
443  }
444}