source: trunk/source/geometry/divisions/test/ExDivisions/src/ExVDivTester.cc @ 1316

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 19.5 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: ExVDivTester.cc,v 1.7 2009/05/16 09:15:33 arce Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30// class ExVDivTester Implementation file
31//
32// 26.05.03 - P.Arce Initial version
33// ********************************************************************
34
35#include "ExVDivTester.hh"
36
37#include "G4Box.hh"
38#include "G4PVPlacement.hh"
39
40#include "Randomize.hh"
41
42#include "G4PVDivisionFactory.hh"
43#include "G4ReflectionFactory.hh"
44#include "G4PVReplica.hh"
45#include "G4PVDivision.hh"
46
47#include "G4Material.hh"
48#include "G4VisAttributes.hh"
49#include "G4Colour.hh"
50
51#include <sstream>
52#include <fstream>
53
54G4bool ExVDivTester::bDivCylindrical = 0;
55
56//--------------------------------------------------------------------------
57ExVDivTester::
58ExVDivTester( PVType& pvtype, PlaceType& postype,
59              std::vector<G4String>& extraPars )
60  : thePVType( pvtype ), thePosType( postype ), theExtraPars( extraPars )
61{
62  verbose = 2;
63  //--- set the number of replicas
64  theWorldLengthXY = 1*m;
65  theWorldLengthZ = 8*theWorldLengthXY;
66  theWorldGap = 10*cm;
67  //--- default material
68  theMate = new G4Material("Pb", 82., 207.19*g/mole, 11.35*g/cm3);
69
70  theDivType = DivNDIV;
71  theOffsetFactor = 0.;
72  theWidthFactor = 1.;
73  theNDiv = 3;
74
75  theStartPhi = 0.*deg;
76  theDeltaPhi = 360.*deg;
77
78  // Create the division factory singleton
79  //
80  G4PVDivisionFactory::GetInstance();
81
82  //---- extract PVType from extraPars (NREPLICAS, WIDTH, NREPLICAS&WIDTH)
83  if( theExtraPars.size() != 0 )
84  {
85    readExtraPars();
86  }
87}
88
89
90//--------------------------------------------------------------------------
91void ExVDivTester::readExtraPars()
92{
93  G4int ii, siz = theExtraPars.size();
94  G4bool ifound;
95  for( ii = 0; ii < siz; ii++ ){
96    ifound = 0;
97    if( theExtraPars[ii].substr(0,6) == "offset" ) {
98      theOffsetFactor = getValueFromExtraPar( theExtraPars[ii] );
99      ifound = 1;
100    }else if( theExtraPars[ii].substr(0,5) == "width" ) {
101      theWidthFactor = getValueFromExtraPar( theExtraPars[ii] );
102      ifound = 1;
103    } else if( theExtraPars[ii].substr(0,7) == "divtype" ) {
104      G4int ival = int(getValueFromExtraPar( theExtraPars[ii] ));
105      if( ival == 0 ){
106        theDivType = DivNDIVandWIDTH;
107        ifound = 1;
108      }else if( ival == 1 ){
109        theDivType = DivNDIV;
110        ifound = 1;
111      }else if( ival == 2 ){
112        theDivType = DivWIDTH;
113        ifound = 1;
114      }
115    }else if( theExtraPars[ii].substr(0,8) == "startphi" ) {
116      theStartPhi = getValueFromExtraPar( theExtraPars[ii] )*deg;
117      ifound = 1;
118    }else if( theExtraPars[ii].substr(0,8) == "deltaphi" ) {
119      theDeltaPhi = getValueFromExtraPar( theExtraPars[ii] )*deg;
120      ifound = 1;
121    }
122
123    if( !ifound ) {
124      G4String msg = "Extra Parameter not in list! : "+theExtraPars[ii];
125      G4Exception("ExVDivTester::readExtraPars()",
126                  "IllegalConstruct", FatalException, msg.c_str() );
127    }
128  }
129}
130
131
132//---------------------------------------------------------------------------
133G4double ExVDivTester::getValueFromExtraPar( G4String& epar )
134{
135  G4double val = 0.;
136
137  G4int ieq = epar.find("=");
138  if( ieq == -1 ) {
139    G4String msg ="Extra Parameter does not have an '='! :"+epar;
140    G4Exception("ExVDivTester::getValueFromExtraPar()",
141                "IllegalConstruct", FatalException, msg.c_str() );
142
143  } else {
144    G4String eparValue = epar.substr(ieq+1,epar.size() );
145    val = atof( eparValue.c_str() );
146  }
147
148  return val;
149}
150
151
152//--------------------------------------------------------------------------
153ExVDivTester::~ExVDivTester()
154{
155  delete theMate;
156}
157
158//--------------------------------------------------------------------------
159void ExVDivTester::GenerateRandomPoints()
160{
161  std::ofstream fout("points.lis");
162  G4double posX, posY, posZ;
163
164  numberOfPoints = 1000;
165  for( G4int ii = 0; ii < numberOfPoints; ii++ ) {
166    posX = CLHEP::RandFlat::shoot(-theWorldLengthXY, theWorldLengthXY);
167    posY = CLHEP::RandFlat::shoot(-theWorldLengthXY, theWorldLengthXY);
168    posZ = CLHEP::RandFlat::shoot(-theWorldLengthZ, theWorldLengthZ);
169    fout << posX << " " << posY << " " << posZ << G4endl;
170  }
171}
172
173//--------------------------------------------------------------------------
174G4VPhysicalVolume* ExVDivTester::BuildGeometry()
175{
176  // Build the world volume
177  //
178  G4Box *worldSolid =
179    new G4Box ("Big Cube", theWorldLengthXY+theWorldGap,
180                           theWorldLengthXY+theWorldGap,
181                           theWorldLengthZ+theWorldGap);
182
183  // Air
184  G4double a, density;
185  a = 14.01*g/mole;
186  G4Element* elN = new G4Element("Nitrogen", "N", 7., a);
187  a = 16.00*g/mole;
188  G4Element* elO = new G4Element("Oxigen", "O", 8., a);
189  density = 1.29*mg/cm3;
190  G4Material* Air = new G4Material("Air", density, 2);
191  Air->AddElement(elN, .7);
192  Air->AddElement(elO, .3);
193
194  G4LogicalVolume *worldLog=new G4LogicalVolume(worldSolid,Air,
195                                                "World",0,0,0);
196 
197  G4PVPlacement *worldPhys=new G4PVPlacement(0,G4ThreeVector(0,0,0),
198                                             "World",worldLog,
199                                             0,false,0);
200  G4VisAttributes* worldVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
201  worldLog->SetVisAttributes(worldVisAtt);
202
203  if( verbose >= 1 )
204    G4cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@ BuildParentSolids "<< G4endl;
205  BuildParentSolids();
206  if( verbose >= 1 )
207    G4cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@ BuildParentVolumes "<< G4endl;
208  BuildParentVolumes( worldLog );
209
210  if( verbose >= 1 )
211    G4cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@ BuildChildrenSolids "<< G4endl;
212  BuildChildrenSolids();
213  if( verbose >= 1 )
214    G4cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@ SetOffsets "<< G4endl;
215  SetOffsets();
216  if( verbose >= 1 )
217    G4cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@ BuildChildrenVolumes "<< G4endl;
218  BuildChildrenVolumes(); //LV and divisions as PV
219
220
221  return worldPhys;
222}
223
224//--------------------------------------------------------------------------
225void ExVDivTester::BuildParentVolumes( G4LogicalVolume* worldLog )
226{
227  // build three volumes
228  G4LogicalVolume* parent1Log =
229    new G4LogicalVolume(theParentSolids[0],theMate,"parentLog-1",0,0,0);
230  G4LogicalVolume* parent2Log =
231    new G4LogicalVolume(theParentSolids[1],theMate,"parentLog-2",0,0,0);
232  G4LogicalVolume* parent3Log =
233    new G4LogicalVolume(theParentSolids[2],theMate,"parentLog-3",0,0,0);
234
235  G4VisAttributes* parentVisAtt = new G4VisAttributes(G4Colour(1.0,0.0,0.0));
236  parent1Log->SetVisAttributes(parentVisAtt);
237  parent2Log->SetVisAttributes(parentVisAtt);
238  parent3Log->SetVisAttributes(parentVisAtt);
239  if ( thePosType == pvReflected )
240  {
241    G4LogicalVolume* parent1LogRefl
242      = G4ReflectionFactory::Instance()->GetReflectedLV(parent1Log);
243    G4LogicalVolume* parent2LogRefl
244      = G4ReflectionFactory::Instance()->GetReflectedLV(parent2Log);
245    G4LogicalVolume* parent3LogRefl
246      = G4ReflectionFactory::Instance()->GetReflectedLV(parent2Log);
247    if(parent1LogRefl) parent1LogRefl->SetVisAttributes(parentVisAtt);
248    if(parent2LogRefl) parent2LogRefl->SetVisAttributes(parentVisAtt);
249    if(parent3LogRefl) parent3LogRefl->SetVisAttributes(parentVisAtt);
250  }
251
252  theParentLogs.push_back(parent1Log);
253  theParentLogs.push_back(parent2Log);
254  theParentLogs.push_back(parent3Log);
255 
256  //parent physical volumes positions do not depend on SolidType, PVType
257  G4int ii;
258  G4int nParents = theAxis.size();
259  for( ii = 0; ii < nParents; ii++ )
260  {
261    G4String parentstr = "parent-";
262    std::ostringstream os;
263    os << ii;
264    G4String buf = os.str();
265    parentstr += buf;
266    if ( thePosType == pvReflected )
267    {
268      // With reflection
269      G4Translate3D translate(0.,0.,(ii-1)*7*theWorldLengthXY);
270      G4ReflectZ3D  reflect;
271      G4Transform3D transform = translate * reflect; 
272
273      G4PhysicalVolumesPair pvPair
274        = G4ReflectionFactory::Instance()
275          ->Place(transform, parentstr, theParentLogs[ii], worldLog, FALSE, 0); 
276      theParentPhyss.push_back(pvPair.first);
277    }
278    else
279    {
280      // Without reflection
281      theParentPhyss.push_back( new G4PVPlacement( 0,
282                              G4ThreeVector(0.,0.,(ii-1)*7*theWorldLengthXY),
283                              theParentLogs[ii] , parentstr,
284                              worldLog, FALSE, 0 ) );
285    }
286  }
287}
288
289//--------------------------------------------------------------------------
290void ExVDivTester::SetOffsets()
291{
292  G4int ii;
293  G4int nParents = theAxis.size();
294  for( ii = 0; ii < nParents; ii++ )
295  {
296    theOffsets.push_back( theWidths[ii]*theOffsetFactor );
297    if( verbose >= 2 )
298      G4cout << ii << " offsets " << theOffsets[ii] << " "
299             << theWidths[ii] << G4endl;
300  }
301}
302
303
304//--------------------------------------------------------------------------
305void ExVDivTester::BuildChildrenVolumes()
306{
307  G4int ii;
308  G4VisAttributes* childVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,0.0));
309  G4VisAttributes* childVisAtt2 = new G4VisAttributes(G4Colour(0.0,1.0,0.0));
310
311  G4int nParents = theAxis.size();
312  for( ii = 0; ii < nParents; ii++ )
313  {
314    G4String childlogstr = "childLog-";
315    std::ostringstream os;
316    os << ii;
317    G4String buf = os.str();
318    childlogstr += buf;
319    G4LogicalVolume* childLog = new G4LogicalVolume(theChildSolids[ii],
320                                theMate, childlogstr,0,0,0);
321    childLog->SetVisAttributes(childVisAtt);
322    theChildLogs.push_back(childLog);
323  }
324
325  //----- Build divisions
326  if( thePVType == pvDivision )
327  {
328    for( ii = 0; ii < nParents; ii++ )
329    {
330      if( verbose >= 1 )
331        G4cout << " @@@@ Building Child volume " << ii << G4endl;
332      G4String childstr = "child-";
333      std::ostringstream os;
334      os << ii;
335      G4String buf = os.str();
336      childstr += buf;
337
338      if ( thePosType == pvReflected )
339      {
340        // With reflection
341        //
342        G4ReflectionFactory::Instance()->SetVerboseLevel(2);
343        if( theDivType == DivNDIVandWIDTH )
344        {
345          G4ReflectionFactory::Instance()
346            ->Divide(childstr, theChildLogs[ii], theParentLogs[ii],
347                               theAxis[ii],
348                               theNDiv, 
349                               theWidths[ii]*theWidthFactor,
350                               theOffsets[ii] );
351
352          if( verbose >= 1 ) G4cout << "division NDIVandWIDTH " << theNDiv
353                 << " " << theWidths[ii]*theWidthFactor
354                 << " " << theOffsets[ii]*theOffsetFactor << G4endl;
355        }
356        else if( theDivType == DivNDIV )
357        {
358          G4ReflectionFactory::Instance()
359            ->Divide(childstr, theChildLogs[ii], theParentLogs[ii],
360                               theAxis[ii],
361                               theNDiv, 
362                               theOffsets[ii] );
363          if( verbose >= 1 )
364            G4cout << "division NDIV " << theNDiv << " "
365                   << theOffsets[ii] << G4endl;
366        }
367        else if( theDivType == DivWIDTH )
368        {
369          G4ReflectionFactory::Instance()
370            ->Divide(childstr, theChildLogs[ii], theParentLogs[ii],
371                               theAxis[ii],
372                               theWidths[ii]*theWidthFactor,
373                               theOffsets[ii] );
374          if( verbose >= 1 )
375            G4cout << "division WIDTH " << theWidths[ii]*theWidthFactor
376                   << " " << theOffsets[ii]*theOffsetFactor << G4endl;
377        }
378
379        // set vis attributes to reflected volumes
380        G4LogicalVolume* childLogRefl
381          = G4ReflectionFactory::Instance()->GetReflectedLV(theChildLogs[ii]);
382        if (childLogRefl) childLogRefl->SetVisAttributes(childVisAtt2);
383      }
384      else
385      {
386        // Without reflection
387        //
388        if( theDivType == DivNDIVandWIDTH )
389        {
390          new G4PVDivision(childstr, theChildLogs[ii], theParentLogs[ii], 
391                           theAxis[ii],
392                           theNDiv, 
393                           theWidths[ii]*theWidthFactor,
394                           theOffsets[ii] );
395          if( verbose >= 1 ) G4cout << "division NDIVandWIDTH " << theNDiv
396                 << " " << theWidths[ii]*theWidthFactor
397                 << " " << theOffsets[ii]*theOffsetFactor << G4endl;
398        }
399        else if( theDivType == DivNDIV )
400        {
401          new G4PVDivision(childstr, theChildLogs[ii], theParentLogs[ii], 
402                           theAxis[ii],
403                           theNDiv, 
404                           theOffsets[ii] );
405          if( verbose >= 1 )
406            G4cout << "division NDIV " << theNDiv << " "
407                   << theOffsets[ii] << G4endl;
408        }
409        else if( theDivType == DivWIDTH )
410        {
411          new G4PVDivision(childstr, theChildLogs[ii], theParentLogs[ii], 
412                           theAxis[ii],
413                           theWidths[ii]*theWidthFactor, 
414                           theOffsets[ii] );
415          if( verbose >= 1 )
416            G4cout << "division WIDTH " << theWidths[ii]*theWidthFactor
417                   << " " << theOffsets[ii]*theOffsetFactor << G4endl;
418        }
419      }
420    }
421  }
422  else if( thePVType == pvReplica )
423  {
424    G4int nParents = theAxis.size();
425    for( ii = 0; ii < nParents; ii++ )
426    {
427      G4String childstr = "child-";
428      std::ostringstream os;
429      os << ii;
430      G4String buf = os.str();
431      childstr += buf;
432      if ( thePosType == pvReflected )
433      {
434        G4ReflectionFactory::Instance()
435          ->Replicate(childstr, theChildLogs[ii], theParentLogs[ii],
436                                theAxis[ii],
437                                theNDiv, 
438                                theWidths[ii]*theWidthFactor,
439                                theOffsets[ii] );
440
441        // set vis attributes to reflected volumes
442        G4LogicalVolume* childLogRefl
443          = G4ReflectionFactory::Instance()->GetReflectedLV(theChildLogs[ii]);
444        if (childLogRefl) childLogRefl->SetVisAttributes(childVisAtt2);
445      }
446      else
447      {
448        new G4PVReplica(childstr, theChildLogs[ii], theParentLogs[ii], 
449                                  theAxis[ii],
450                                  theNDiv, 
451                                  theWidths[ii]*theWidthFactor, 
452                                  theOffsets[ii] );
453      }
454      if( verbose >= 1 )
455        G4cout << "replica " <<  theNDiv << " "
456               << theWidths[ii]*theWidthFactor << " "
457               << theOffsets[ii]*theOffsetFactor << G4endl;
458    }
459  }
460}
461
462//------------------------------------------------------------------------
463void ExVDivTester::PrintChildrenSolids( std::ostream& out )
464{
465  for( size_t ii = 0; ii < theChildSolids.size(); ii++)
466  {
467    out << theChildSolids[0] << "child solid built  0 "
468        << *theChildSolids[0] << G4endl;
469    out << " theChildSolid built " <<  ii << G4endl;
470    out <<  *theChildSolids[ii] << G4endl;
471  }
472}
473
474//------------------------------------------------------------------------
475void ExVDivTester::PrintParentSolid( std::ostream& out )
476{
477  out << " PARENT SOLID " << G4endl;
478
479  G4int nParents = theAxis.size();
480  for( G4int ii = 0; ii < nParents; ii++ ){
481    out << *(theParentSolids[ii]) << G4endl;
482  }
483}
484
485
486//--------------------------------------------------------------------------
487void ExVDivTester::GenerateScanPointsAsBox()
488{
489  std::ofstream fout("points.lis");
490  G4int ii;
491
492  G4int nPointsPerDiv = 2;
493  G4double disppoint = 0*mm;
494  numberOfPoints = theNDiv * nPointsPerDiv;
495  // For division along X
496  G4ThreeVector centre(0.,0.,-theWorldLengthZ+theWorldLengthXY);
497  for( ii = 0; ii < numberOfPoints; ii++ )
498  {
499    // any Z, any Y
500    G4ThreeVector pR( 0., disppoint, disppoint );
501    G4double X = -theWorldLengthXY
502               + (ii+0.5) * 2*theWorldLengthXY/numberOfPoints;
503    pR.setX( X );
504    pR += centre;
505    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
506  }
507
508  // For division along Y
509  centre = G4ThreeVector(0.,0.,0.);
510  for( ii = 0; ii < numberOfPoints; ii++ )
511  {
512    // any X, any Z
513    G4ThreeVector pR(disppoint, 0., disppoint );
514    G4double Y = -theWorldLengthXY
515               + (ii+0.5) * 2*theWorldLengthXY/numberOfPoints;
516    pR.setY( Y );
517    pR += centre;
518    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
519  }
520
521
522  // For division along Z
523  centre = G4ThreeVector(0.,0.,theWorldLengthZ-theWorldLengthXY);
524  for( ii = 0; ii < numberOfPoints; ii++ )
525  {
526    // any X, any Y
527    G4ThreeVector pR( disppoint, disppoint, 0. );
528    G4double Z = -theWorldLengthXY
529               + (ii+0.5) * 2*theWorldLengthXY/numberOfPoints;
530    pR.setZ( Z );
531    pR += centre;
532    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
533  }
534}
535
536
537//--------------------------------------------------------------------------
538void ExVDivTester::GenerateScanPointsAsTubs()
539{
540  std::ofstream fout("points.lis");
541  G4int ii;
542
543  G4int nPointsPerDiv = 5;
544  numberOfPoints = theNDiv*nPointsPerDiv;
545  if( verbose >= 1 ) G4cout << " numberOfPoints " << numberOfPoints << G4endl;
546  // For division along R
547  G4ThreeVector centre(0.,0.,-theWorldLengthZ+theWorldLengthXY);
548  for( ii = 0; ii < numberOfPoints; ii++ )
549  {
550    // any Z, phi close initial phi
551    G4double phi = 0.*deg;
552    //  if( theExtraPars.size() > 0 ) phi = (theExtraPars[0]+1.)*deg;
553    G4ThreeVector pR( std::cos(phi), std::sin(phi), 0. );
554    G4double rho = (ii+0.5) * theWorldLengthXY/numberOfPoints;
555    pR.setRho( rho );
556    pR += centre;
557    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
558  }
559 
560  // For division along phi
561  centre = G4ThreeVector(theWorldLengthXY*0.7,0.,0.);
562  for( ii = 0; ii < numberOfPoints; ii++ )
563  {
564    G4double phi = (ii+0.5) * 360*deg/numberOfPoints;
565    // any Z
566    //    G4ThreeVector pR( std::cos(phi), std::sin(phi), theWorldLengthXY/100. );
567    //    pR += centre;
568    G4ThreeVector pR( centre );
569    pR.setPhi( phi );
570    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
571  }
572 
573  // For division along Z
574  centre = G4ThreeVector(theWorldLengthXY*0.7,
575                         0., theWorldLengthZ-theWorldLengthXY);
576  for( ii = 0; ii < numberOfPoints; ii++ )
577  {
578    //any rho (=1), phi close initial phi
579    G4ThreeVector pR( centre);
580    G4double Z = centre.z() - theWorldLengthXY
581               + (ii+0.5) * 2*theWorldLengthXY / numberOfPoints;
582    pR.setZ( Z );
583    fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl;
584  }
585}
Note: See TracBrowser for help on using the repository browser.