DetectorConstruction.cc 17.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// Author: Michael Phipps

#include "DetectorConstruction.hh"
28
#include "DetectorMessenger.hh"
29
30
31
32
33
34
35

#include "G4GeometryManager.hh"
#include "G4SolidStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4PhysicalVolumeStore.hh"

#include "G4RunManager.hh"
36
#include "G4UImanager.hh"
37
38
39
40
41
42
43
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"

44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "G4EqMagElectricField.hh"
#include "PurgMagTabulatedField3D.hh"
#include "G4ChordFinder.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4PVParameterised.hh"
#include "G4PVReplica.hh"
#include "G4UniformMagField.hh"
#include "G4ExplicitEuler.hh"
#include "G4ImplicitEuler.hh"
#include "G4SimpleRunge.hh"
#include "G4SimpleHeum.hh"
#include "G4ClassicalRK4.hh"
#include "G4HelixExplicitEuler.hh"
#include "G4HelixImplicitEuler.hh"
#include "G4HelixSimpleRunge.hh"
#include "G4CashKarpRKF45.hh"
#include "G4RKG3_Stepper.hh"


65
66
67
68
69
70
#include <iostream>
#include <stdio.h>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

DetectorConstruction::DetectorConstruction()
71
  : G4VUserDetectorConstruction(), m_solidWorld(NULL), m_logicWorld(NULL),
72
  m_physWorld(NULL),CLUSTER(false),OPTICAL(false), ForceDetectorPosition(false)
73
{
74
  new DetectorMessenger(this);
75
76
77
78
  currentRPD = -1;
  currentZDC = -1;
  m_materials = Materials::getInstance();
}
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

DetectorConstruction::~DetectorConstruction()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* DetectorConstruction::Construct(){
  if ( m_physWorld ) {
    G4GeometryManager::GetInstance()->OpenGeometry();
    G4PhysicalVolumeStore::GetInstance()->Clean();
    G4LogicalVolumeStore::GetInstance()->Clean();
    G4SolidStore::GetInstance()->Clean();
    G4LogicalSkinSurface::CleanSurfaceTable();
    G4LogicalBorderSurface::CleanSurfaceTable();
  }

97
98
99
  G4UImanager* UImanager = G4UImanager::GetUIpointer();
  UImanager->ApplyCommand("/control/execute geometry.mac");

100
101
102
  m_materials->UseOpticalMaterials(OPTICAL);
  m_materials->DefineOpticalProperties();

103
104
105
106
  if( ForceDetectorPosition ){
    ManualConstruction();
  }else{
    ConstructSPSTestBeam();
107
    std::cout << "built test beam" << std::endl;
108
109
  }
  return m_physWorld;
110
}
111

112
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113

114
115
G4VPhysicalVolume* DetectorConstruction::ConstructWorldVolume(G4double x, G4double y, G4double z){
  printf( "Building world with x %5.1f y %5.1f z %5.1f\n", x, y, z );
116
117

  m_solidWorld =
118
    new G4Box("World", 0.5*x, 0.5*y, 0.5*z );
119
120
121

  m_logicWorld =
    new G4LogicalVolume(m_solidWorld,     //its solid
122
                        m_materials->Air, //its material
123
                        "World");         //its name
124
  std::cout << "The world's name is " << m_logicWorld->GetName() <<  std::endl;
125
126
127
128
129
130
131
132
133
134

  m_physWorld =
    new G4PVPlacement(0,                  //no rotation
                      G4ThreeVector(),    //at (0,0,0)
                      m_logicWorld,       //its logical volume
                      "World",            //its name
                      0,                  //its mother  volume
                      false,              //no boolean operation
                      0,                  //copy number
                      false);             //overlaps checking
135

136

137
  G4VisAttributes* boxVisAtt_world = new G4VisAttributes(G4Colour(0.5,0.5,0.5));
138

139
140
  m_logicWorld ->SetVisAttributes(boxVisAtt_world);

141
142
143
144
145
146
147
148
149
150
151
152
  return m_physWorld;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* DetectorConstruction::ManualConstruction(){

  std::cout << "******************************************" << std::endl
            << "        PLACING DETECTORS MANUALLY        " << std::endl
            << "******************************************" << std::endl;

153
154
155
156
157
158
159
160
161
162
163
164
  G4ThreeVector* pos;
  for(ModTypeZDC* zdc : m_ZDCvec){
    pos = zdc->GetPosition();
    printf( "ZDC%d center = (%f,%f,%f)\n", zdc->GetModNum(), pos->x(), pos->y(), pos->z() );
    zdc->Construct();
  }

  for(ModTypeRPD* rpd : m_RPDvec){
    pos = rpd->GetPosition();
    printf( "RPD%d center = (%f,%f,%f)", rpd->GetModNum(), pos->x(), pos->y(), pos->z() );
    rpd->Construct();
  }
165

166
  return m_physWorld;
167
}
168
169


170
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171

172
173
174
175
176
177
178
179
180
181
182
183
G4VPhysicalVolume* DetectorConstruction::ConstructSPSTestBeam(){
  // Create variables to be used in beamtest 2018 simulation
  G4ThreeVector zdc1Pos,zdc2Pos, rpdPos;
  bool ZDC1 = false, ZDC2 = false, RPD = false;
  G4double mag_zOffset=0, firstDetZ, detX, detY, detZ,
           tableX_shift=0, tableY_shift=0;
  G4Material* Lead = m_materials->Pb;
  G4double worldSizeX = 180*mm;
  G4double worldSizeY = 1200*mm;
  G4double worldSizeZ = 32000*mm;
  G4double worldZoffset= worldSizeZ/2.0;
  std::string detector[3];
184

185
  //################################ World volume construction
186

187
  ConstructWorldVolume( worldSizeX, worldSizeY, worldSizeZ );
188

189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  //################################ SURVEY/ALIGNMENT_SETUP

  LoadConfigurationFile();
  LoadAlignmentFile();

  //table(-2250,500) -> rpd/beam(0,0)	where 100=0.1cm in table coordinates
  //-320mm is offset to get from zdc mid to active area mid
  tableX_shift = (-2250.0 - (m_alignment->x_table)  )/100*mm ;//2257 more accurate
  tableY_shift = (500.0   - (m_alignment->y_table)  )/100*mm ;//501  more accurate
  Survey *srvy_rpd = GetSurvey("RPD");

  detector[0]=m_alignment->upstream_Det;
  detector[1]=m_alignment->mid_Det;
  detector[2]=m_alignment->downstream_Det;

  for(int i=0; i<3; i++){
    if(detector[i]=="ZDC1") {ZDC1 = true;}
    if(detector[i]=="ZDC2") {ZDC2 = true;}
    if(detector[i]=="RPD" ) {RPD  = true;}
  }

  firstDetZ = worldSizeZ;
  for(Survey* survey : m_surveyEntries){
    detX = (survey->x_pos*1000.0)*mm;
    detY = (survey->y_pos*1000.0)*mm;
    detZ = (survey->z_pos*1000.0 - worldZoffset )*mm;

    if(detZ < firstDetZ) firstDetZ = detZ;

    if( survey->detector == "ZDC1" && ZDC1 ){
      AddZDC( new G4ThreeVector( detX + (      - srvy_rpd->x_pos + tableX_shift)*mm,
                                 detY + (- 230 - srvy_rpd->y_pos + tableY_shift)*mm,
                                 detZ ) );
    } else if( survey->detector == "ZDC2" && ZDC2 ){
      AddZDC( new G4ThreeVector( detX + (      - srvy_rpd->x_pos + tableX_shift)*mm,
                                 detY + (- 230 - srvy_rpd->y_pos + tableY_shift)*mm,
                                 detZ ) );
    }else if( survey->detector == "RPD" && RPD ){
      AddRPD( new G4ThreeVector( detX, detY, detZ) );
    }
229
  }
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245

  G4ThreeVector* pos;
  G4int modNum;
  for(ModTypeZDC* zdc : m_ZDCvec){
    zdc->SetFiberDiameters    ( new G4ThreeVector(1.5*mm,0.0,0.0) );
    zdc->SetAbsorberDimensions( new G4ThreeVector(90.0*mm, 180.0*mm, 11.0*mm) );
    zdc->SetnAbsorbers        ( 11 );
    zdc->SetHousingThickness  ( 4.5*mm );
    zdc->SetGapThickness      ( 2.5*mm );
    zdc->SetHousingMaterial   ( "aluminum" );
    zdc->SetAbsorberMaterial  ( "pure" );

    pos = zdc->GetPosition();
    modNum = zdc->GetModNum();
    printf( "ZDC%d center = (%f,%f,%f)", modNum, pos->x(), pos->y(), pos->z() );
    zdc->Construct();
246
  }
247
248

  for(ModTypeRPD* rpd : m_RPDvec){
249
250
    rpd->SetFiberDiameters( new G4ThreeVector(0.6*mm,0.68*mm,0.72*mm) );
    rpd->SetDetectorType  ( "cms" );
251
252
253
254
255

    pos = rpd->GetPosition();
    modNum = rpd->GetModNum();
    printf( "RPD%d center = (%f,%f,%f)", modNum, pos->x(), pos->y(), pos->z() );
    rpd->Construct();
256
  }
257
//################################ Get SURVEY/ALIGNMENT_END
258

259

260
  // Setup magnetic field
261
  if( m_alignment->magnet_On ){
262
263
264
265
    mag_zOffset=-9.55*m;
    //Field grid in A9.TABLE. File must be accessible from run directory.
    G4MagneticField* PurgMagField= new PurgMagTabulatedField3D((std::getenv("JZCaPA") + std::string("/Utils/PurgMag3D.TABLE")).c_str(), mag_zOffset+(worldZoffset/1000.0));
    fField.Put(PurgMagField);
266

267
268
    //This is thread-local
    G4FieldManager* pFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
269

270
271
272
273
    pFieldMgr->SetDetectorField(fField.Get());
    pFieldMgr->CreateChordFinder(fField.Get());
  }
    // Setup lead target
274
    if( m_alignment->target_In ){
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
      G4Box* leadTarget = new G4Box("Target",5*cm, 5*cm, 1.3*cm);

      logic_leadTarget
      = new G4LogicalVolume(leadTarget,
                            Lead,
                            "LeadTarget");

      new G4PVPlacement(0,
                        G4ThreeVector(0.0, 0.0, (2600-worldZoffset)*mm),
                        logic_leadTarget,
                        "LeadTarget1",
                        m_logicWorld,
                        false,
                        0
                        );

    // Visualization attributes
    G4VisAttributes* boxVisAtt_lead= new G4VisAttributes(G4Colour(1.0,0.0,1.0)); //magenta
    logic_leadTarget ->SetVisAttributes(boxVisAtt_lead);
294

295
   }
296
297
298
299
300
301
   if( m_alignment->lead_In ){
     // Assign lead block position in mm, place approx 1 ft in front of ZDC1
     // half ZDC Z width = 90mm
     // 1ft ~ 300mm
     G4double leadblockZ = ( firstDetZ - 90 - 300)*mm;

302
     G4Box* leadBlock = new G4Box("LeadBlock",(0.5*worldSizeX)*mm, (0.5*150)*mm, 10*cm);
303

304
305
306
307
     logic_leadBlock
     = new G4LogicalVolume(leadBlock,
                           Lead,
                           "LeadBlock");
308

309
    new G4PVPlacement(0,
aricct2's avatar
aricct2 committed
310
                      G4ThreeVector(0.0, 0.0, (leadblockZ)*mm),
311
312
313
314
315
316
317
                      logic_leadBlock,
                      "LeadBlock1",
                      m_logicWorld,
                      false,
                      0
                      );

318
319
320
321
322
    // Visualization attributes
    G4VisAttributes* boxVisAtt_leadblk= new G4VisAttributes(G4Colour(1.0,0.0,1.0)); //magenta
    logic_leadBlock ->SetVisAttributes(boxVisAtt_leadblk);
    cout << "Placed Lead Block at Z = " << leadblockZ*mm << "mm" <<  std::endl;
  }
323
324
  return m_physWorld;
}
325
326
327
328
329
330
331
332
333
334
335
336
337





//***********************************************************************************
// READ Survey_2018.xml and Alignment_2018.xml
//***********************************************************************************
/**
 * @brief Reads the .xml configuration file and load characteristics for all the alignments, immediately sorted into alignment objects
 * @param _inFile
 */

338
339
340
341
void DetectorConstruction::LoadConfigurationFile( G4String _inFile  ){

  if(_inFile = ""){
    _inFile = std::getenv("JZCaPA");
342
    _inFile += "/Utils/Survey_2018.xml";
343
  }
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378

    m_XMLparser = new XMLSettingsReader();


    if (!m_XMLparser->parseFile(_inFile)) {
      if(!m_XMLparser->parseFile("Survey_2018.xml")){
        std::cerr << " Data Reader could not parse file : " << _inFile << std::endl;
        return;
      }
      std::cout << " Found Survey_2018.xml in executable directory " << std::endl;
    }

    std::cout << "Loading .xml Configuration File..." << std::endl;
    std::cout << "Found " << m_XMLparser->getBaseNodeCount("Survey") << " survey entries " << std::endl;

    int first_run, last_run;

    for (int i = 0; i < m_XMLparser->getBaseNodeCount("Survey"); i++) { //this was unsigned int i = 0

        m_XMLparser->getChildValue("Survey",i,"start_run",first_run);
        m_XMLparser->getChildValue("Survey",i,"end_run",last_run);

        //Discard entries for any alignment that does not apply to our run
        if(m_runNumber < first_run || m_runNumber > last_run) continue;
          m_survey = new Survey();

        //If the entry applies, we store it in the vector
        m_XMLparser->getChildValue("Survey",i,"detector",m_survey->detector);
        m_XMLparser->getChildValue("Survey",i,"x_pos",m_survey->x_pos);
        m_XMLparser->getChildValue("Survey",i,"y_pos",m_survey->y_pos);
        m_XMLparser->getChildValue("Survey",i,"z_pos",m_survey->z_pos);
        m_XMLparser->getChildValue("Survey",i,"cos_x",m_survey->cos_x);
        m_XMLparser->getChildValue("Survey",i,"cos_y",m_survey->cos_y);
        m_XMLparser->getChildValue("Survey",i,"cos_z",m_survey->cos_z);

379
        m_surveyEntries.push_back(m_survey);
380
381
    }

382
    if(m_surveyEntries.size() == 0) std::cout << "WARNING: SURVEY NOT FOUND!!!" << std::endl;
383
384
385
386
387
388
389
390

    return;
}

/**
 * @brief Reads the .xml configuration file and load characteristics for all the channels, immediately sorted into detectors objects
 * @param _inFile
 */
391
void DetectorConstruction::LoadAlignmentFile( G4String _inFile ){
392
393
  bool debug = false;

394
395
  if( _inFile = ""){
    _inFile = std::getenv("JZCaPA");
396
    _inFile += "/Utils/Alignment_2018.xml";
397
398
  }

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
    m_XMLparser = new XMLSettingsReader();

    if (!m_XMLparser->parseFile(_inFile)) {
        if(!m_XMLparser->parseFile("Alignment_2018.xml")){
          std::cerr << " Data Reader could not parse file : " << _inFile << std::endl;
          return;
        }
        std::cout << " Found Alignment_2018.xml in executable directory " << std::endl;
    }

    m_alignment = new Alignment();

    if(debug){
      std::cout << "Loading .xml Alignment File..." << std::endl;
      std::cout << "Found " << m_XMLparser->getBaseNodeCount("Alignment") << " alignment entries " << std::endl;
      std::cout << "Retrieving the information for run " << m_runNumber << std::endl;
    }

    int run;
    for ( int i = 0; i < m_XMLparser->getBaseNodeCount("Alignment"); i++) {
      m_XMLparser->getChildValue("Alignment",i,"run",run);
      if(run != m_runNumber) continue;
        if(debug){
          std::cout << "Found Run Entry in Alignment file for run " << m_runNumber << std::endl;
        }
        m_XMLparser->getChildValue("Alignment",i,"x_table",m_alignment->x_table);
        m_XMLparser->getChildValue("Alignment",i,"y_table",m_alignment->y_table);
        m_XMLparser->getChildValue("Alignment",i,"upstream_Det",m_alignment->upstream_Det);
        m_XMLparser->getChildValue("Alignment",i,"mid_Det",m_alignment->mid_Det);
        m_XMLparser->getChildValue("Alignment",i,"downstream_Det",m_alignment->downstream_Det);
        m_XMLparser->getChildValue("Alignment",i,"target_In",m_alignment->target_In);
        m_XMLparser->getChildValue("Alignment",i,"lead_In",m_alignment->lead_In);
        m_XMLparser->getChildValue("Alignment",i,"magnet_On",m_alignment->magnet_On);
    }

    if(m_alignment == NULL) std::cout << "WARNING: ALIGNMENT NOT FOUND!!!" << std::endl;
    return;
}

438
439
440
441
442
443
/*
 * Retrieve a survey based on name
*/
Survey* DetectorConstruction::GetSurvey(G4String name){
  for(unsigned int i = 0; i < m_surveyEntries.size(); i++){
    if( m_surveyEntries[i]->detector.compare( name.c_str() ) ){ return m_surveyEntries[i]; }
444
445
446
447
448
  }
  Survey* empty=NULL;
  return empty;
}

449
450
451
452
/*
 * Add a ZDC module
*/
void DetectorConstruction::AddZDC(G4ThreeVector* position){
453
  uint newModNum = m_ZDCvec.size()+1;
454
  m_ZDCvec.push_back(new ModTypeZDC(newModNum, m_logicWorld, position ));
455
456
  currentZDC = newModNum;
  printf("Added ZDC%d\n", newModNum);
457
458
459
460
461
462
}

/*
 * Add an RPD module
*/
void DetectorConstruction::AddRPD(G4ThreeVector* position){
463
  uint newModNum = m_RPDvec.size()+1;
464
  m_RPDvec.push_back(new ModTypeRPD(newModNum, m_logicWorld, position ));
465
466
  currentRPD = newModNum;
  printf("Added RPD%d\n", newModNum);
467
468
469
470
471
472
}

/*
 * Duplicate a ZDC module
*/
void DetectorConstruction::DuplicateZDC( G4int module ){
473
474
  uint newModNum = m_ZDCvec.size()+1;
  if((unsigned)module >= newModNum ){
475
476
    printf("\n\n Cannot duplicate. ZDC%d does not exist \n\n",module);
  }
477
478
479
  m_ZDCvec.push_back( new ModTypeZDC( newModNum, m_ZDCvec.at(module-1) ) );
  currentZDC = newModNum;
  printf("Duplicate ZDC%d built from ZDC%d\n", newModNum, module);
480
481
482
483
484
485
}

/*
 * Duplicate an RPD module
*/
void DetectorConstruction::DuplicateRPD( G4int module ){
486
487
  uint newModNum = m_RPDvec.size()+1;
  if((unsigned)module >= newModNum ){
488
489
    printf("\n\n Cannot duplicate. RPD%d does not exist \n\n",module);
  }
490
491
492
  m_RPDvec.push_back( new ModTypeRPD( newModNum, m_RPDvec.at(module-1) ) );
  currentRPD = newModNum;
  printf("Duplicate RPD%d built from RPD%d\n", newModNum, module);
493
}