source: trunk/source/geometry/divisions/src/G4ParameterisationTrd.cc

Last change on this file was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 16.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4ParameterisationTrd.cc,v 1.19 2010/11/10 09:16:08 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// class G4ParameterisationTrd Implementation file
31//
32// 26.05.03 - P.Arce, Initial version
33// 08.04.04 - I.Hrivnacova, Implemented reflection
34// 21.04.10 - M.Asai, Added gaps
35// --------------------------------------------------------------------
36
37#include "G4ParameterisationTrd.hh"
38
39#include <iomanip>
40#include "G4ThreeVector.hh"
41#include "G4RotationMatrix.hh"
42#include "G4VPhysicalVolume.hh"
43#include "G4LogicalVolume.hh"
44#include "G4ReflectedSolid.hh"
45#include "G4Trd.hh"
46#include "G4Trap.hh"
47
48//--------------------------------------------------------------------------
49G4VParameterisationTrd::
50G4VParameterisationTrd( EAxis axis, G4int nDiv, G4double width,
51                        G4double offset, G4VSolid* msolid,
52                        DivisionType divType )
53  :  G4VDivisionParameterisation( axis, nDiv, width, offset, divType, msolid ),
54     bDivInTrap(false)
55{
56  G4Trd* msol = (G4Trd*)(msolid);
57  if (msolid->GetEntityType() == "G4ReflectedSolid")
58  {
59    // Get constituent solid 
60    G4VSolid* mConstituentSolid
61       = ((G4ReflectedSolid*)msolid)->GetConstituentMovedSolid();
62    msol = (G4Trd*)(mConstituentSolid);
63 
64    // Create a new solid with inversed parameters
65    G4Trd* newSolid
66      = new G4Trd(msol->GetName(),
67                  msol->GetXHalfLength2(), msol->GetXHalfLength1(),
68                  msol->GetYHalfLength2(), msol->GetYHalfLength1(),
69                  msol->GetZHalfLength());
70    msol = newSolid;
71    fmotherSolid = newSolid;
72    fReflectedSolid = true;
73    fDeleteSolid = true;
74  }   
75}
76
77//------------------------------------------------------------------------
78G4VParameterisationTrd::~G4VParameterisationTrd()
79{
80}
81
82//------------------------------------------------------------------------
83G4ParameterisationTrdX::
84G4ParameterisationTrdX( EAxis axis, G4int nDiv,
85                        G4double width, G4double offset,
86                        G4VSolid* msolid, DivisionType divType )
87  :  G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType )
88{
89  CheckParametersValidity();
90  SetType( "DivisionTrdX" );
91
92  G4Trd* msol = (G4Trd*)(fmotherSolid);
93  if( divType == DivWIDTH )
94  {
95    fnDiv = CalculateNDiv( msol->GetXHalfLength1()+msol->GetXHalfLength2(),
96                           width, offset );
97  }
98  else if( divType == DivNDIV )
99  {
100    fwidth = CalculateWidth( msol->GetXHalfLength1()+msol->GetXHalfLength2(), nDiv, offset );
101  }
102
103#ifdef G4DIVDEBUG
104  if( verbose >= 1 )
105  {
106    G4cout << " G4ParameterisationTrdX - ## divisions " << fnDiv << " = "
107           << nDiv << G4endl
108           << " Offset " << foffset << " = " << offset << G4endl
109           << " Width " << fwidth << " = " << width << G4endl;
110  }
111#endif
112
113  G4double mpDx1 = msol->GetXHalfLength1();
114  G4double mpDx2 = msol->GetXHalfLength2();
115  if( std::fabs(mpDx1 - mpDx2) > kCarTolerance )
116  {
117    bDivInTrap = true;
118  }
119}
120
121//------------------------------------------------------------------------
122G4ParameterisationTrdX::~G4ParameterisationTrdX()
123{
124}
125
126//------------------------------------------------------------------------
127G4double G4ParameterisationTrdX::GetMaxParameter() const
128{
129  G4Trd* msol = (G4Trd*)(fmotherSolid);
130  return (msol->GetXHalfLength1()+msol->GetXHalfLength2());
131}
132
133//------------------------------------------------------------------------
134void
135G4ParameterisationTrdX::
136ComputeTransformation( const G4int copyNo,
137                       G4VPhysicalVolume *physVol ) const
138{
139  G4Trd* msol = (G4Trd*)(fmotherSolid );
140  G4double mdx = ( msol->GetXHalfLength1() + msol->GetXHalfLength2() ) / 2.;
141
142  //----- translation
143  G4ThreeVector origin(0.,0.,0.); 
144  G4double posi;
145  if( !bDivInTrap )
146  {
147    posi = -mdx + foffset + (copyNo+0.5)*fwidth;
148  }
149  else
150  {
151    G4double aveHL = (msol->GetXHalfLength1()+msol->GetXHalfLength2())/2.;
152    posi = - aveHL + foffset + (copyNo+0.5)*aveHL/fnDiv*2;
153  }
154  if( faxis == kXAxis )
155  {
156    origin.setX( posi ); 
157  }
158  else
159  { 
160    G4Exception("G4ParameterisationTrdX::ComputeTransformation()",
161                "IllegalConstruct", FatalException,
162                "Only axes along X are allowed, and axis is: "+faxis);
163  }
164
165#ifdef G4DIVDEBUG
166  if( verbose >= 2 )
167  {
168    G4cout << std::setprecision(8) << " G4ParameterisationTrdX::ComputeTransformation "
169           << copyNo << G4endl
170           << " Position: " << origin << " - Axis: " << faxis << G4endl;
171  }
172#endif
173
174   //----- set translation
175  physVol->SetTranslation( origin );
176}
177
178//--------------------------------------------------------------------------
179void
180G4ParameterisationTrdX::
181ComputeDimensions( G4Trd& trd, const G4int, const G4VPhysicalVolume* ) const
182{
183  G4Trd* msol = (G4Trd*)(fmotherSolid);
184
185  G4double pDy1 = msol->GetYHalfLength1();
186  G4double pDy2 = msol->GetYHalfLength2();
187  G4double pDz = msol->GetZHalfLength();
188  G4double pDx = fwidth/2. - fhgap;
189 
190  trd.SetAllParameters ( pDx, pDx, pDy1, pDy2, pDz );
191
192#ifdef G4DIVDEBUG
193  if( verbose >= 2 )
194  {
195    G4cout << " G4ParameterisationTrdX::ComputeDimensions():" << copyNo << G4endl;
196     trd.DumpInfo();
197  }
198#endif
199}
200
201G4VSolid* G4ParameterisationTrdX::
202ComputeSolid(const G4int i, G4VPhysicalVolume * pv)
203{
204  if( bDivInTrap ) 
205  {
206    return G4VDivisionParameterisation::ComputeSolid(i, pv);
207  } 
208  else 
209  {
210    return fmotherSolid;
211  }
212}
213
214
215//--------------------------------------------------------------------------
216void
217G4ParameterisationTrdX::
218ComputeDimensions( G4Trap& trap, const G4int copyNo, const G4VPhysicalVolume* ) const
219{
220  G4Trd* msol = (G4Trd*)(fmotherSolid);
221  G4double pDy1 = msol->GetYHalfLength1();
222  G4double pDy2 = msol->GetYHalfLength2();
223  G4double pDz = msol->GetZHalfLength();
224  G4double pDx1 = msol->GetXHalfLength1()/fnDiv; 
225  G4double pDx2 = msol->GetXHalfLength2()/fnDiv; 
226
227  G4double cxy1 = -msol->GetXHalfLength1() + foffset + (copyNo+0.5)*pDx1*2;// centre of the side at y=-pDy1
228  G4double cxy2 = -msol->GetXHalfLength2() + foffset + (copyNo+0.5)*pDx2*2;// centre of the side at y=+pDy1
229  G4double alp = std::atan( (cxy2-cxy1)/pDz );
230 
231  trap.SetAllParameters ( pDz,
232                          0.,
233                          0.,
234                          pDy1,
235                          pDx1,
236                          pDx2,
237                          alp,
238                          pDy2,
239                          pDx1 - fhgap,
240                          pDx2 - fhgap * pDx2/pDx1,
241                          alp);
242
243#ifdef G4DIVDEBUG
244  if( verbose >= 2 )
245  {
246    G4cout << " G4ParameterisationTrdX::ComputeDimensions():" << copyNo << G4endl;
247     trap.DumpInfo();
248  }
249#endif
250}
251
252//--------------------------------------------------------------------------
253void G4ParameterisationTrdX::CheckParametersValidity()
254{
255  G4VDivisionParameterisation::CheckParametersValidity();
256/*
257  G4Trd* msol = (G4Trd*)(fmotherSolid);
258
259  G4double mpDx1 = msol->GetXHalfLength1();
260  G4double mpDx2 = msol->GetXHalfLength2();
261  bDivInTrap = false;
262
263  if( std::fabs(mpDx1 - mpDx2) > kCarTolerance )
264  {
265    G4cerr << "ERROR - G4ParameterisationTrdX::CheckParametersValidity()"
266           << G4endl
267           << "        Making a division of a TRD along axis X," << G4endl
268           << "        while the X half lengths are not equal," << G4endl
269           << "        is not (yet) supported. It will result" << G4endl
270           << "        in non-equal division solids." << G4endl;
271    G4Exception("G4ParameterisationTrdX::CheckParametersValidity()",
272                "IllegalConstruct", FatalException,
273                "Invalid solid specification. NOT supported.");
274  }
275*/
276}
277
278//--------------------------------------------------------------------------
279void G4ParameterisationTrdX::
280ComputeTrapParams()
281{
282}
283
284//--------------------------------------------------------------------------
285G4ParameterisationTrdY::
286G4ParameterisationTrdY( EAxis axis, G4int nDiv,
287                        G4double width, G4double offset,
288                        G4VSolid* msolid, DivisionType divType)
289  : G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType )
290{
291  CheckParametersValidity();
292  SetType( "DivisionTrdY" );
293
294  G4Trd* msol = (G4Trd*)(fmotherSolid);
295  if( divType == DivWIDTH )
296  {
297    fnDiv = CalculateNDiv( 2*msol->GetYHalfLength1(),
298                           width, offset );
299  }
300  else if( divType == DivNDIV )
301  {
302    fwidth = CalculateWidth( 2*msol->GetYHalfLength1(),
303                             nDiv, offset );
304  }
305
306#ifdef G4DIVDEBUG
307  if( verbose >= 1 )
308  {
309     G4cout << " G4ParameterisationTrdY no divisions " << fnDiv
310            << " = " << nDiv << G4endl
311            << " Offset " << foffset << " = " << offset << G4endl
312            << " width " << fwidth << " = " << width << G4endl;
313  }
314#endif
315}
316
317//------------------------------------------------------------------------
318G4ParameterisationTrdY::~G4ParameterisationTrdY()
319{
320}
321
322//------------------------------------------------------------------------
323G4double G4ParameterisationTrdY::GetMaxParameter() const
324{
325  G4Trd* msol = (G4Trd*)(fmotherSolid);
326  return 2*msol->GetYHalfLength1();
327}
328
329//--------------------------------------------------------------------------
330void
331G4ParameterisationTrdY::
332ComputeTransformation( const G4int copyNo, G4VPhysicalVolume* physVol ) const
333{
334  G4Trd* msol = (G4Trd*)(fmotherSolid );
335  G4double mdy = msol->GetYHalfLength1();
336
337  //----- translation
338  G4ThreeVector origin(0.,0.,0.); 
339  G4double posi = -mdy + foffset + (copyNo+0.5)*fwidth;
340
341  if( faxis == kYAxis )
342  {
343    origin.setY( posi ); 
344  }
345  else
346  { 
347    G4Exception("G4ParameterisationTrdY::ComputeTransformation()",
348                "IllegalConstruct", FatalException,
349                "Only axes along Y are allowed !");
350  }
351
352#ifdef G4DIVDEBUG
353  if( verbose >= 2 )
354  {
355    G4cout << std::setprecision(8) << " G4ParameterisationTrdY::ComputeTransformation " << copyNo
356           << " pos " << origin << " rot mat " << " axis " << faxis << G4endl;
357  }
358#endif
359
360   //----- set translation
361  physVol->SetTranslation( origin );
362}
363
364//--------------------------------------------------------------------------
365void
366G4ParameterisationTrdY::
367ComputeDimensions(G4Trd& trd, const G4int, const G4VPhysicalVolume*) const
368{
369  //---- The division along Y of a Trd will result a Trd, only
370  //--- if Y at -Z and +Z are equal, else use the G4Trap version
371  G4Trd* msol = (G4Trd*)(fmotherSolid);
372 
373  G4double pDx1 = msol->GetXHalfLength1();
374  G4double pDx2 = msol->GetXHalfLength2();
375  G4double pDz = msol->GetZHalfLength();
376  G4double pDy = fwidth/2. - fhgap;
377 
378  trd.SetAllParameters ( pDx1, pDx2, pDy, pDy, pDz );
379
380#ifdef G4DIVDEBUG
381  if( verbose >= 2 )
382  {
383    G4cout << " G4ParameterisationTrdY::ComputeDimensions():" << G4endl;
384    trd.DumpInfo();
385  }
386#endif
387}
388
389//--------------------------------------------------------------------------
390void G4ParameterisationTrdY::CheckParametersValidity()
391{
392  G4VDivisionParameterisation::CheckParametersValidity();
393
394  G4Trd* msol = (G4Trd*)(fmotherSolid);
395
396  G4double mpDy1 = msol->GetYHalfLength1();
397  G4double mpDy2 = msol->GetYHalfLength2();
398
399  if( std::fabs(mpDy1 - mpDy2) > kCarTolerance )
400  {
401    G4cerr << "ERROR - G4ParameterisationTrdY::CheckParametersValidity()"
402           << G4endl
403           << "        Making a division of a TRD along axis Y while" << G4endl
404           << "        the Y half lengths are not equal is not (yet)" << G4endl
405           << "        supported. It will result in non-equal" << G4endl
406           << "        division solids." << G4endl;
407    G4Exception("G4ParameterisationTrdY::CheckParametersValidity()",
408                "IllegalConstruct", FatalException,
409                "Invalid solid specification. NOT supported.");
410  }
411}
412
413//--------------------------------------------------------------------------
414G4ParameterisationTrdZ::
415G4ParameterisationTrdZ( EAxis axis, G4int nDiv,
416                        G4double width, G4double offset,
417                        G4VSolid* msolid, DivisionType divType )
418  : G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType )
419{ 
420  CheckParametersValidity();
421  SetType( "DivTrdZ" );
422
423  G4Trd* msol = (G4Trd*)(fmotherSolid);
424  if( divType == DivWIDTH )
425  {
426    fnDiv = CalculateNDiv( 2*msol->GetZHalfLength(),
427                           width, offset );
428  }
429  else if( divType == DivNDIV )
430  {
431    fwidth = CalculateWidth( 2*msol->GetZHalfLength(),
432                             nDiv, offset );
433  }
434
435#ifdef G4DIVDEBUG
436  if( verbose >= 1 )
437  {
438    G4cout << " G4ParameterisationTrdZ no divisions " << fnDiv
439           << " = " << nDiv << G4endl
440           << " Offset " << foffset << " = " << offset << G4endl
441           << " Width " << fwidth << " = " << width << G4endl;
442  }
443#endif
444}
445
446//------------------------------------------------------------------------
447G4ParameterisationTrdZ::~G4ParameterisationTrdZ()
448{
449}
450
451//------------------------------------------------------------------------
452G4double G4ParameterisationTrdZ::GetMaxParameter() const
453{
454  G4Trd* msol = (G4Trd*)(fmotherSolid);
455  return 2*msol->GetZHalfLength();
456}
457
458//--------------------------------------------------------------------------
459void
460G4ParameterisationTrdZ::
461ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol) const
462{
463  G4Trd* msol = (G4Trd*)(fmotherSolid );
464  G4double mdz = msol->GetZHalfLength();
465
466  //----- translation
467  G4ThreeVector origin(0.,0.,0.); 
468  G4double posi = -mdz + OffsetZ() + (copyNo+0.5)*fwidth;
469  if( faxis == kZAxis )
470  {
471    origin.setZ( posi ); 
472  }
473  else
474  { 
475    G4Exception("G4ParameterisationTrdZ::ComputeTransformation()",
476                "IllegalConstruct", FatalException,
477                "Only axes along Z are allowed !");
478  }
479
480#ifdef G4DIVDEBUG
481  if( verbose >= 1 )
482  {
483    G4cout << std::setprecision(8) << " G4ParameterisationTrdZ: "
484           << copyNo << G4endl
485           << " Position: " << origin << " - Offset: " << foffset
486           << " - Width: " << fwidth << " Axis " << faxis << G4endl;
487  }
488#endif
489
490   //----- set translation
491  physVol->SetTranslation( origin );
492}
493
494//--------------------------------------------------------------------------
495void
496G4ParameterisationTrdZ::
497ComputeDimensions(G4Trd& trd, const G4int copyNo,
498                  const G4VPhysicalVolume*) const
499{
500  //---- The division along Z of a Trd will result a Trd
501  G4Trd* msol = (G4Trd*)(fmotherSolid);
502
503  G4double pDx1 = msol->GetXHalfLength1();
504  G4double DDx = (msol->GetXHalfLength2() - msol->GetXHalfLength1() );
505  G4double pDy1 = msol->GetYHalfLength1();
506  G4double DDy = (msol->GetYHalfLength2() - msol->GetYHalfLength1() );
507  G4double pDz = fwidth/2. - fhgap;
508  G4double zLength = 2*msol->GetZHalfLength();
509 
510  trd.SetAllParameters( pDx1+DDx*(OffsetZ()+copyNo*fwidth+fhgap)/zLength,
511                        pDx1+DDx*(OffsetZ()+(copyNo+1)*fwidth-fhgap)/zLength, 
512                        pDy1+DDy*(OffsetZ()+copyNo*fwidth+fhgap)/zLength,
513                        pDy1+DDy*(OffsetZ()+(copyNo+1)*fwidth-fhgap)/zLength, pDz );
514
515#ifdef G4DIVDEBUG
516  if( verbose >= 1 )
517  {
518    G4cout << " G4ParameterisationTrdZ::ComputeDimensions()"
519           << " - Mother TRD " << G4endl;
520    msol->DumpInfo();
521    G4cout << " - Parameterised TRD: "
522           << copyNo << G4endl;
523    trd.DumpInfo();
524  }
525#endif
526}
Note: See TracBrowser for help on using the repository browser.