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

Last change on this file since 1202 was 1010, checked in by garnier, 15 years ago

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