source: trunk/source/geometry/magneticfield/test/NTST/src/NTSTDetectorConstruction.cc@ 1306

Last change on this file since 1306 was 1231, checked in by garnier, 16 years ago

update from CVS
File size: 16.6 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#include "NTSTDetectorConstruction.hh"
28#include "NTSTDetectorMessenger.hh"
29#include "NTSTRotationMatrix.hh"
30#include "G4TransportationManager.hh"
31#include "G4FieldManager.hh"
32#include "G4ChordFinder.hh"
33#include "G4Material.hh"
34#include "G4Box.hh"
35#include "G4Tubs.hh"
36#include "G4Trd.hh"
37#include "G4LogicalVolume.hh"
38#include "G4ThreeVector.hh"
39#include "G4PVPlacement.hh"
40#include "G4VisAttributes.hh"
41#include "G4Color.hh"
42#include "G4Transform3D.hh"
43#include "G4Point3D.hh"
44#include "globals.hh"
45#include "NTSTFileRead.hh"
46
47#include <iomanip>
48
49#include "G4Mag_UsualEqRhs.hh"
50#include "G4ClassicalRK4.hh"
51#include "G4SimpleRunge.hh"
52#include "G4CashKarpRKF45.hh"
53#include "G4RKG3_Stepper.hh"
54#include "G4HelixMixedStepper.hh"
55#include "G4NystromRK4.hh"
56
57#include "G4DELPHIMagField.hh"
58#include "G4PropagatorInField.hh"
59
60NTSTDetectorConstruction::NTSTDetectorConstruction()
61 : _FileRead(0), debug(false), radius(19*cm), NSubLayer(0),
62 disableSVT(false), disableDCH(false),
63 field( 1.5*tesla, 0, 0 ),
64 fpChordFinder( 0 ),
65 fMinChordStep( 0.1 ) // was 0.001 *mm )
66{
67 _FileRead = new NTSTFileRead("SVT.dat");
68
69 // create commands necessary for the definition of the SVT
70
71 DetectorMessenger = new NTSTDetectorMessenger(this);
72}
73
74NTSTDetectorConstruction::~NTSTDetectorConstruction()
75{
76 delete _FileRead;
77 delete fpChordFinder;
78 delete DetectorMessenger;
79}
80
81void NTSTDetectorConstruction::SetInputFileName(G4String FileName)
82{
83 delete _FileRead;
84 _FileRead = new NTSTFileRead(FileName);
85}
86
87
88void NTSTDetectorConstruction::SetDebugCmd(G4int NewDebug)
89{
90 debug = NewDebug;
91 if (debug) {
92 G4cout << "Reset debug flag to true" << G4endl;
93 }
94 else {
95 G4cout << "Reset debug flag to false" << G4endl;
96 }
97}
98
99
100void
101NTSTDetectorConstruction::SetNSubLayer(G4int NewNSubLayer)
102{
103 NSubLayer = NewNSubLayer;
104 G4cout << "Reset number of sublayers to " << NSubLayer << G4endl;
105}
106
107void
108NTSTDetectorConstruction::SetOuterRadius(G4double NewRadius)
109{
110 radius = NewRadius;
111 G4cout << "Reset SVT mother volume outer radius parameter to " << radius
112 << G4endl;
113}
114
115void
116NTSTDetectorConstruction::DisableDetector(G4String theDetector)
117{
118 if (theDetector == "SVT") {
119 G4cout << "Disable " << theDetector << " detector" << G4endl;
120 disableSVT=true;
121 }
122 else if (theDetector == "DCH") {
123 G4cout << "Disable " << theDetector << " detector" << G4endl;
124 disableDCH=true;
125 }
126 else if (theDetector == "all") {
127 G4cout << "Disable SVT and DCH" << G4endl;
128 disableSVT=true;
129 disableDCH=true;
130 }
131 else if (theDetector == "none") {
132 G4cout << "Enable SVT and DCH" << G4endl;
133 disableSVT=false;
134 disableDCH=false;
135 }
136}
137
138
139void
140NTSTDetectorConstruction::PrintCorners(const G4Transform3D& theT,
141 G4LogicalVolume* theLV)
142{
143 G4VSolid* theSolid=theLV->GetSolid();
144 G4Trd* theTRD=(G4Trd*)theSolid;
145
146 G4double x1=theTRD->GetXHalfLength1();
147 G4double x2=theTRD->GetXHalfLength2();
148 G4double y1=theTRD->GetYHalfLength1();
149 G4double y2=theTRD->GetYHalfLength2();
150 G4double z =theTRD->GetZHalfLength();
151
152 G4Point3D t1(-x1, -y1, -z);
153 G4Point3D t2(+x1, -y1, -z);
154 G4Point3D t3(-x1, +y1, -z);
155 G4Point3D t4(+x1, +y1, -z);
156 G4Point3D t5(-x2, -y2, z);
157 G4Point3D t6(+x2, -y2, z);
158 G4Point3D t7(-x2, +y2, z);
159 G4Point3D t8(+x2, +y2, z);
160
161 G4Point3D u1 = theT*t1;
162 G4Point3D u2 = theT*t2;
163 G4Point3D u3 = theT*t3;
164 G4Point3D u4 = theT*t4;
165 G4Point3D u5 = theT*t5;
166 G4Point3D u6 = theT*t6;
167 G4Point3D u7 = theT*t7;
168 G4Point3D u8 = theT*t8;
169
170 G4cout << std::setw(9) << u1.z() << std::setw(9) << u2.z() << std::setw(9) << u3.z()
171 << std::setw(9) << u4.z() << std::setw(9) << u5.z() << std::setw(9) << u6.z()
172 << std::setw(9) << u7.z() << std::setw(9) << u8.z() << G4endl;
173}
174
175
176G4VPhysicalVolume*
177NTSTDetectorConstruction::Construct()
178{
179 //------------------------------------------------------ field
180 G4Mag_UsualEqRhs *pEquation;
181 G4MagIntegratorStepper *pStepper;
182 G4FieldManager *globalFieldManager;
183
184 globalFieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager();
185
186 G4PropagatorInField *
187 globalPropagatorInField= G4TransportationManager::GetTransportationManager()->GetPropagatorInField();
188
189 globalPropagatorInField->SetMaxLoopCount( 10000 );
190 G4cout
191 << "PropagatorInField parameter(s) are: " << G4endl
192 << " SetMaxLoopCount=" << globalPropagatorInField->GetMaxLoopCount()
193 << " minEpsilonStep= " << globalPropagatorInField->GetMinimumEpsilonStep() << " "
194 << " maxEpsilonStep= " << globalPropagatorInField->GetMaximumEpsilonStep() << " "
195 << G4endl;
196
197 globalFieldManager->SetDetectorField( (G4MagneticField *)&field );
198
199 // globalFieldManager->SetMinimumEpsilonStep( 5.0e-7 ); // Old value
200 // globalFieldManager->SetMaximumEpsilonStep( 0.05 ); // FIX - old value
201 // globalFieldManager->SetDeltaOneStep( 0.25 * mm ); // original value
202 // globalFieldManager->SetDeltaIntersection( 0.10 * mm ); // original value
203
204 G4cout << "Field Manager's parameters are "
205 << " minEpsilonStep= " << globalFieldManager->GetMinimumEpsilonStep() << " "
206 << " maxEpsilonStep= " << globalFieldManager->GetMaximumEpsilonStep() << " "
207 << " deltaOneStep= " << globalFieldManager->GetDeltaOneStep() << " "
208 << " deltaIntersection= " << globalFieldManager->GetDeltaIntersection()
209 << G4endl;
210
211 pEquation = new G4Mag_UsualEqRhs( &field);
212
213 // pStepper =
214 // new G4ClassicalRK4( pEquation ); G4cout << "Stepper is " << "ClassicalRK4" << G4endl;
215 // new G4RKG3_Stepper( pEquation ); // Nystrom, like Geant3
216 // pStepper= new G4SimpleRunge( pEquation ); G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
217 // pStepper= new G4CashKarpRKF45( pEquation ); G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
218 // pStepper= new G4HelixMixedStepper( pEquation ); G4cout << "Stepper is " << "HelixMixed" << G4endl;
219 // pStepper= StepperFactory::CreateStepper( order );
220
221 pStepper= new G4NystromRK4( pEquation ); G4cout << "Stepper is " << "NystromRK4" << G4endl;
222
223 // G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
224 // << "ClassicalRK4" << G4endl;
225 // << " G4HelixMixedStepper " << G4endl;
226
227 // globalFieldManager->CreateChordFinder( (G4MagneticField *)&field );
228 fpChordFinder= new G4ChordFinder( (G4MagneticField *)&field,
229 fMinChordStep,
230 pStepper );
231 fpChordFinder->SetVerbose(1);
232 globalFieldManager->SetChordFinder( fpChordFinder );
233
234 //------------------------------------------------------ materials
235
236 G4double a; // atomic mass
237 G4double z; // atomic number
238 G4double density;
239 G4String name;
240
241 a = 39.95*g/mole;
242 density = 1.782e-03*g/cm3;
243 G4Material* Ar = new G4Material(name="ArgonGas", z=18., a, density);
244
245 a = 26.98*g/mole;
246 density = 2.7*g/cm3;
247 // G4Material* Al =
248 new G4Material(name="Aluminum", z=13., a, density);
249
250 a = 28.0855*g/mole;
251 density = 2.33*g/cm3;
252 G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
253
254 //------------------------------------------------------ volumes
255
256 //------------------------------ experimental hall (world volume)
257
258 G4double expHall_x = 1000*cm;
259 G4double expHall_y = 1000*cm;
260 G4double expHall_z = 2000*cm;
261 G4Box* experimentalHall_box
262 = new G4Box("expHall_box",expHall_x,expHall_y,expHall_z);
263 G4LogicalVolume* experimentalHall_log
264 = new G4LogicalVolume(experimentalHall_box,Ar,"expHall_log",0,0,0);
265
266 experimentalHall_log->SetVisAttributes(G4VisAttributes::Invisible);
267
268 G4VPhysicalVolume* experimentalHall_phys
269 = new G4PVPlacement(0,G4ThreeVector(),"expHall",
270 experimentalHall_log,0,false,0);
271
272 G4double innerRadiusOfTheSvt = 2.9*cm;
273 G4double outerRadiusOfTheSvt = radius;
274 G4double lengthOfTheSvt = 40.*cm;
275 G4double startAngleOfTheSvt = 0*deg;
276 G4double spanningAngleOfTheSvt = 360.*deg;
277
278 G4double SvtPos_x = 0.*m;
279 G4double SvtPos_y = 0.*m;
280 G4double SvtPos_z = 0.*m;
281
282 G4double innerRadiusOfTheDch = 24*cm;
283 G4double outerRadiusOfTheDch = 81*cm;
284 G4double lengthOfTheDch = 250*cm;
285 G4double startAngleOfTheDch = 0*deg;
286 G4double spanningAngleOfTheDch = 360.*deg;
287
288 G4double DchPos_x = 0.*m;
289 G4double DchPos_y = 0.*m;
290 G4double DchPos_z = 0.*m;
291
292 disableSVT=false;
293
294 if (disableSVT == false){
295
296 //------------------------------ SVT tracker volume
297
298 G4Tubs* Svt_tube
299 = new G4Tubs("Svt_tube",innerRadiusOfTheSvt,
300 outerRadiusOfTheSvt,lengthOfTheSvt,
301 startAngleOfTheSvt,spanningAngleOfTheSvt);
302 G4LogicalVolume* Svt_log
303 = new G4LogicalVolume(Svt_tube,Ar,"Svt_log",0,0,0);
304
305 Svt_log -> SetVisAttributes(G4VisAttributes::Invisible);
306
307 // G4VPhysicalVolume* Svt_phys =
308 new G4PVPlacement(0,
309 G4ThreeVector(SvtPos_x,
310 SvtPos_y,
311 SvtPos_z),
312 Svt_log,"Svt",experimentalHall_log,false,0);
313
314 if (debug)
315 G4cout << "Placed SVT mother of length: "
316 << std::setw(7) << lengthOfTheSvt/cm
317 << " and radii (cm): "
318 << std::setw(7) << innerRadiusOfTheSvt/cm
319 << std::setw(7) << outerRadiusOfTheSvt/cm << G4endl;
320
321 //------------------------------ SVT guts
322
323 // read in parameters of the wafers
324
325 int NwafType=0;
326 _FileRead->StreamLine() >> NwafType;
327 if (debug) G4cout << "Number of wafer types: " << NwafType << G4endl;
328
329 G4LogicalVolume** theWafer_log = new G4LogicalVolume*[NwafType];
330
331 // define wafer vis attributes
332
333 G4Color red(1,0,0);
334 // G4Color green(0,1,0);
335 // G4Color blue(0,0,1);
336 G4VisAttributes* vAttr = new G4VisAttributes(red);
337 // make solid
338 vAttr->SetForceSolid(true);
339
340 // define wafer shapes and create logical volumes indexed by Wafer type
341
342 for (int ind=0; ind<NwafType; ind++){
343
344 G4double Hmin, Hmax, Hzlen, Hthick;
345 G4int IwafType;
346 _FileRead->StreamLine() >> IwafType >> Hmin >> Hmax >> Hzlen >> Hthick;
347 if (debug) G4cout << "Wafer type " << std::setw(3) << IwafType
348 << " Hmin " << std::setw(10) << Hmin/cm
349 << " Hmax " << std::setw(10) << Hmax/cm
350 << " Hzlen " << std::setw(10) << Hzlen/cm
351 << " Hthick " << std::setw(6) << Hthick/cm << G4endl;
352
353 G4Trd* aWafer = new G4Trd("aWafer", Hthick*mm, Hthick*mm,
354 Hmin*mm, Hmax*mm, Hzlen*mm);
355 theWafer_log[IwafType-1]
356 = new G4LogicalVolume(aWafer, Si, "aWafer_log", 0,0,0);
357
358 theWafer_log[IwafType-1] -> SetVisAttributes(vAttr);
359 }
360
361 // get number of layers
362
363 G4int Nsublayer=NSubLayer;
364
365 _FileRead->StreamLine() >> Nsublayer;
366 if (debug) G4cout << "Number of layers " << Nsublayer << G4endl;
367
368 if (NSubLayer>0 && NSubLayer<=7){
369 Nsublayer = NSubLayer;
370 }
371
372 // loop over the number of layers
373
374 for (G4int Isublay=0; Isublay<Nsublayer;Isublay++){
375 G4int Ilayer, Isublayer, Nmodule;
376 _FileRead->StreamLine() >> Ilayer >> Isublayer >> Nmodule;
377 if (debug) G4cout << "Number of modules for layer "
378 << std::setw(3) << Ilayer
379 << " sublayer " << std::setw(3) << Isublayer << " = "
380 << std::setw(3) << Nmodule << G4endl;
381
382 // loop over the number of modules
383
384 for (G4int Imod=0; Imod<Nmodule; Imod++){
385 G4int Imodule, Nwafer;
386 _FileRead->StreamLine() >> Imodule >> Nwafer;
387 if (debug) G4cout << "Number of wafers in module "
388 << std::setw(3) << Imodule
389 << " = " << std::setw(3) << Nwafer << G4endl;
390
391 // loop over the number of wafers in a module
392
393 for (G4int Iwaf=0; Iwaf < Nwafer; Iwaf++){
394 G4int Iwafer, IwaferType;
395 _FileRead->StreamLine() >> Iwafer >> IwaferType;
396 if (debug) G4cout << "Wafer " << std::setw(3) << Iwafer
397 << " type " << std::setw(3)
398 << IwaferType << G4endl;
399 G4double x,y,z;
400 _FileRead->StreamLine() >> x >> y >> z;
401 G4ThreeVector WafPos(x*mm, y*mm, z*mm);
402 if (debug) G4cout << " position " << std::setw(9) << x << " "
403 << std::setw(9) << y
404 << " " << std::setw(9) << z << G4endl;
405
406 _FileRead->StreamLine() >> x >> y >> z;
407 if (debug) G4cout << "Rotation Matrix:" << G4endl;
408
409 G4ThreeVector row1(x,y,z);
410 if (debug) G4cout << row1 << G4endl;
411
412 _FileRead->StreamLine() >> x >> y >> z;
413
414 G4ThreeVector row2(x,y,z);
415 if (debug) G4cout << row2 << G4endl;
416
417 _FileRead->StreamLine() >> x >> y >> z;
418
419 G4ThreeVector row3(x,y,z);
420 if (debug) G4cout << row3 << G4endl;
421
422 NTSTRotationMatrix WafMat;
423
424 WafMat.SetRotationMatrixByRow(row1,row2,row3);
425
426 G4Transform3D theTransform(WafMat, WafPos);
427
428 // G4VPhysicalVolume* wafer_phys =
429 new G4PVPlacement(theTransform, theWafer_log[IwaferType-1],
430 "WaferPos",Svt_log,false,0);
431 if (Imod==0 && debug) {
432 G4cout << "lay " << std::setw(3) << Ilayer << " Waf "
433 << Iwafer;
434 PrintCorners(theTransform, theWafer_log[IwaferType-1]);
435 }
436 }
437 }
438 }
439 } // end SVT block
440 if (disableDCH == false) {
441 G4Tubs* Dch_tube
442 = new G4Tubs("Dch_tube",innerRadiusOfTheDch,
443 outerRadiusOfTheDch,lengthOfTheDch,
444 startAngleOfTheDch,spanningAngleOfTheDch);
445 G4LogicalVolume* Dch_log
446 = new G4LogicalVolume(Dch_tube,Ar,"Dch_log",0,0,0);
447
448 Dch_log -> SetVisAttributes(G4VisAttributes::Invisible);
449
450 // G4VPhysicalVolume* Dch_phys =
451 new G4PVPlacement(0,
452 G4ThreeVector(DchPos_x,
453 DchPos_y,
454 DchPos_z),
455 Dch_log,"Dch",experimentalHall_log,false,0);
456
457 if (debug)
458 G4cout << "Placed DCH mother of length: "
459 << std::setw(7) << lengthOfTheDch/cm
460 << " and radii (cm): "
461 << std::setw(7) << innerRadiusOfTheDch/cm
462 << std::setw(7) << outerRadiusOfTheDch/cm << G4endl;
463
464 G4double r[41] = {25, 26, 27, 28, 30, 32, 33, 34, 35, 37, 38, 39, 41, 42,
465 43, 45, 46, 48, 49, 50, 52, 53, 54, 56, 57, 59, 60, 61,
466 62, 64, 66, 67, 68, 70, 71, 72, 73, 75, 76, 77, 78
467 };
468
469 for (int lay=0; lay < 40; lay++){
470 G4double innerRadiusOfTheLayer=r[lay]*cm;
471 G4double outerRadiusOfTheLayer=r[lay+1]*cm;
472 G4double lengthOfTheLayer=lengthOfTheDch;
473 G4double startAngleOfTheLayer=0*deg;
474 G4double spanningAngleOfTheLayer=360*deg;
475
476 G4Tubs* LayTub
477 = new G4Tubs("Lay_tube",innerRadiusOfTheLayer,
478 outerRadiusOfTheLayer,lengthOfTheLayer,
479 startAngleOfTheLayer,spanningAngleOfTheLayer);
480 G4LogicalVolume* Layer_log
481 = new G4LogicalVolume(LayTub,Ar,"Layer_log",0,0,0);
482
483 // G4VPhysicalVolume* Layer_phys =
484 new G4PVPlacement(0,
485 G4ThreeVector(0),
486 Layer_log,"Layer", Dch_log,false,0);
487 if (debug)
488 G4cout << "Placed LAYER mother of length: "
489 << std::setw(7) << lengthOfTheLayer/cm
490 << " and radii (cm): "
491 << std::setw(7) << innerRadiusOfTheLayer/cm
492 << std::setw(7) << outerRadiusOfTheLayer/cm << G4endl;
493 }
494
495 } // end DCH block
496
497
498 //------------------------------------------------------------------
499
500 return experimentalHall_phys;
501}
502
503
Note: See TracBrowser for help on using the repository browser.