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

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