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

Last change on this file since 1207 was 1199, checked in by garnier, 16 years ago

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