Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/divisions/src/G4ParameterisationTrd.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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