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

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

geant4 tag 9.4
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-ref-00 $
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.