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

Last change on this file since 1348 was 1231, checked in by garnier, 15 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.