EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ProjCrystalCalorimeterDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ProjCrystalCalorimeterDetector.cc
1 #include "PHG4ProjCrystalCalorimeterDetector.h"
2 
3 #include "PHG4CrystalCalorimeterDetector.h"
4 #include "PHG4CrystalCalorimeterDisplayAction.h"
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <phool/recoConsts.h>
9 
10 #include <Geant4/G4Cons.hh>
11 #include <Geant4/G4GenericTrap.hh>
12 #include <Geant4/G4LogicalVolume.hh>
13 #include <Geant4/G4Material.hh>
14 #include <Geant4/G4PVPlacement.hh>
15 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
16 #include <Geant4/G4String.hh> // for G4String
17 #include <Geant4/G4SubtractionSolid.hh>
18 #include <Geant4/G4SystemOfUnits.hh>
19 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
20 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
21 #include <Geant4/G4Trd.hh>
22 #include <Geant4/G4TwoVector.hh>
23 #include <Geant4/G4Types.hh> // for G4double, G4int
24 
25 #include <TSystem.h>
26 
27 #include <cmath> // for M_PI
28 #include <cstdlib>
29 #include <iostream>
30 #include <sstream>
31 #include <vector> // for vector
32 
33 class G4VSolid;
34 class PHCompositeNode;
35 
36 using namespace std;
37 
38 //_______________________________________________________________________
39 PHG4ProjCrystalCalorimeterDetector::PHG4ProjCrystalCalorimeterDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
40  : PHG4CrystalCalorimeterDetector(subsys, Node, parameters, dnam)
41  // _dx_front(50.19*mm), //****************************************************************//
42  // _dy_front(50.19*mm), //****************************************************************//
43  // _dx_back(59.3154545455*mm), // PANDA eEMCAL Numbers: Crystals are 2.4cm * 2.4cm on front face //
44  // _dy_back(59.3154545455*mm), //****************************************************************//
45  // _dz_crystal(90.000*mm), //****************************************************************//
46  , _dx_front(41.44 * mm)
47  , _dy_front(41.44 * mm)
48  , _dx_back(48.97454545455 * mm)
49  , _dy_back(48.97454545455 * mm)
50  , _dz_crystal(90.000 * mm)
51  , _crystallogicnameprefix("eEcalCrystal")
52  , m_IsActive(GetParams()->get_int_param("active"))
53  , m_AbsorberActive(GetParams()->get_int_param("absorberactive"))
54 {
55 }
56 
57 //_______________________________________________________________________
58 int PHG4ProjCrystalCalorimeterDetector::IsInCrystalCalorimeter(G4VPhysicalVolume *volume) const
59 {
60  // if hit is in absorber material
61  // bool isinabsorber = false;
62 
63  if (m_IsActive)
64  {
65  if (m_ActiveVolumeSet.find(volume) != m_ActiveVolumeSet.end())
66  {
67  return GetCaloType();
68  }
69  }
70  if (m_AbsorberActive)
71  {
72  if (m_PassiveVolumeSet.find(volume) != m_PassiveVolumeSet.end())
73  {
74  return -1;
75  }
76  }
77  return 0;
78 }
79 
80 //_______________________________________________________________________
81 void PHG4ProjCrystalCalorimeterDetector::ConstructMe(G4LogicalVolume *logicWorld)
82 {
83  if (Verbosity() > 0)
84  {
85  cout << "PHG4ProjCrystalCalorimeterDetector: Begin Construction" << endl;
86  }
87 
88  if (GetParams()->get_string_param("mappingtower").empty())
89  {
90  cout << "ERROR in PHG4ProjCrystalCalorimeterDetector: No tower mapping file specified. Abort detector construction." << endl;
91  cout << "Please run set_string_param(\"mappingtower\", std::string filename ) first." << endl;
92  exit(1);
93  }
94 
95  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
96  recoConsts *rc = recoConsts::instance();
97  G4Material *WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
98 
99  G4VSolid *ecal_envelope_cone = new G4Cons("eEcal_envelope_solid",
100  GetParams()->get_double_param("rMin1") * cm, GetParams()->get_double_param("rMax1") * cm,
101  GetParams()->get_double_param("rMin2") * cm, GetParams()->get_double_param("rMax2") * cm,
102  GetParams()->get_double_param("dz") * cm / 2.,
103  0, 2 * M_PI);
104 
105  G4LogicalVolume *ecal_envelope_log = new G4LogicalVolume(ecal_envelope_cone, WorldMaterial, G4String("eEcal_envelope"), 0, 0, 0);
106  GetDisplayAction()->AddVolume(ecal_envelope_log, "Envelope");
107  /* Define rotation attributes for envelope cone */
108  G4RotationMatrix ecal_rotm;
109  ecal_rotm.rotateX(GetParams()->get_double_param("rot_x") * deg);
110  ecal_rotm.rotateY(GetParams()->get_double_param("rot_y") * deg);
111  ecal_rotm.rotateZ(GetParams()->get_double_param("rot_z") * deg);
112 
113  /* Place envelope cone in simulation */
114  new G4PVPlacement(G4Transform3D(ecal_rotm, G4ThreeVector(GetParams()->get_double_param("place_x") * cm, GetParams()->get_double_param("place_y") * cm, GetParams()->get_double_param("place_z") * cm)),
115  ecal_envelope_log, "CrystalCalorimeter", logicWorld, 0, false, OverlapCheck());
116 
117  /* Construct crystal calorimeter within envelope */
118  ConstructProjectiveCrystals(ecal_envelope_log);
119 
120  return;
121 }
122 
123 void PHG4ProjCrystalCalorimeterDetector::GetCarbonFiberAdjustments(G4double &adjust_width, G4double &adjust_length)
124 {
125  adjust_width = 0.1258525627 * mm; //Because the crystals are slightly angled, the carbon fiber needs to be shortened
126  adjust_length = 2.4824474402 * mm; // from the mother volume (to prevent clipping) by this amount.
127 }
128 
129 void PHG4ProjCrystalCalorimeterDetector::GetCarbonFiberSpacing(G4double &CF_width, G4double &Air_CF, G4double &Air_Cry)
130 {
131  //Parameters of the spacing given by PANDA document arXiv:0810.1216v1 Fig. 7.25
132  CF_width = 0.18 * mm; //Width of the carbon fiber which surrounds the crystal
133  Air_CF = 0.24 * mm; //Air gap between crystal and the carbon fiber
134  Air_Cry = 0.60 * mm; //Air gap between crystal and crystal
135 }
136 
137 int PHG4ProjCrystalCalorimeterDetector::Fill4x4Unit(G4LogicalVolume *crystal_logic)
138 {
139  //*************************************
140  //**********Define Materials***********
141  //*************************************
142 
143  //Crystal Material (Default is Lead Tungstate)
144  G4Material *material_crystal = GetDetectorMaterial(GetParams()->get_string_param("material"));
145 
146  recoConsts *rc = recoConsts::instance();
147  G4Material *WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
148 
149  //*************************************
150  //**********Build First Crystal********
151  //*************************************
152 
153  //Crystal Dimensions determined by the _dx_front, with various gaps and carbon fiber widths subtracted out
154 
155  //Parameters of the spacing given by PANDA document arXiv:0810.1216v1 Fig. 7.25
156  G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
157  GetCarbonFiberSpacing(carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals);
158 
159  //Crystal Dimensions
160  G4double dx_front_small = (_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the front crystal face
161  G4double dy_front_small = (_dy_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full height of the front crystal face
162  G4double dx_back_small = (_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the back crystal face
163  G4double dy_back_small = (_dy_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full height of the back crystal face
164  G4double dz = _dz_crystal;
165 
166  //Vertices of the primary, irregularly shaped crystal put into a vector
167  std::vector<G4TwoVector> vertices;
168  vertices.push_back(G4TwoVector(0, 0));
169  vertices.push_back(G4TwoVector(0, dy_front_small));
170  vertices.push_back(G4TwoVector(dx_front_small, dy_front_small));
171  vertices.push_back(G4TwoVector(dx_front_small, 0));
172  vertices.push_back(G4TwoVector(0, 0));
173  vertices.push_back(G4TwoVector(0, dy_back_small));
174  vertices.push_back(G4TwoVector(dx_back_small, dy_back_small));
175  vertices.push_back(G4TwoVector(dx_back_small, 0));
176 
177  //Create the primary, irregularly shaped crystal
178  G4VSolid *crystal_solid_small = new G4GenericTrap(G4String("eEcal_crystal"),
179  dz,
180  vertices);
181 
182  G4LogicalVolume *crystal_logic_small = new G4LogicalVolume(crystal_solid_small,
183  material_crystal,
184  "eEcal_crystal",
185  0, 0, 0);
186 
187  GetDisplayAction()->AddVolume(crystal_logic_small, "Crystal");
188 
189  //****************************************************
190  //Build the solid/logical volume for the 2 x 2 crystal
191  //****************************************************
192 
193  G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
194  G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
195  G4double TwoByTwo_dy1 = TwoByTwo_dx1;
196  G4double TwoByTwo_dy2 = TwoByTwo_dx2;
197  G4double TwoByTwo_dz = _dz_crystal;
198 
199  G4VSolid *Two_by_Two_solid = new G4Trd(G4String("Two_by_Two_solid"),
200  TwoByTwo_dx1, //Half length on the small face in x
201  TwoByTwo_dx2, //Half length on the large face in x
202  TwoByTwo_dy1, //Half length on the small face in y
203  TwoByTwo_dy2, //Half length on the large face in y
204  TwoByTwo_dz); //Half length in z
205 
206  G4LogicalVolume *Two_by_Two_logic = new G4LogicalVolume(Two_by_Two_solid,
207  WorldMaterial,
208  "2_by_2_unit",
209  0, 0, 0);
210 
211  GetDisplayAction()->AddVolume(Two_by_Two_logic, "TwoByTwo");
212 
213  //*************************************************************
214  //Read in mapping file for a single 2 x 2 block and 4 x 4 block
215  //*************************************************************
216 
217  //The first four lines of the data file refer to the 2x2 block, and the last four lines refer to the mapping of the 4x4 block
218 
219  const int NumberOfIndices = 9; //Number of indices in mapping file for 4x4 block
220 
221  //Find the number of lines in the file, make and fill a NumberOfLines by NumberOfIndices matrix with contents of data file
222  int NumberOfLines = 0;
223  ifstream in(GetParams()->get_string_param("mapping4x4"));
224  if (!in.is_open())
225  {
226  cout << endl
227  << "*******************************************************************" << endl;
228  cout << "ERROR in 2 by 2 crystal mapping ";
229  cout << "Failed to open " << GetParams()->get_string_param("mapping4x4") << " --- Exiting program." << endl;
230  cout << "*******************************************************************" << endl
231  << endl;
232  gSystem->Exit(1);
233  }
234  std::string unused;
235  while (std::getline(in, unused))
236  {
237  ++NumberOfLines;
238  }
239  in.close();
240 
241  in.open(GetParams()->get_string_param("mapping4x4"));
242  G4int j_cry = NumberOfLines;
243  G4int k_cry = NumberOfIndices;
244 
245  double TwoByTwo[j_cry][k_cry];
246 
247  G4int j = 0;
248  G4int k = 0;
249 
250  while (j_cry > j)
251  {
252  while (k_cry > k)
253  {
254  in >> TwoByTwo[j][k];
255  k++;
256  }
257  j++;
258  k = 0;
259  }
260  in.close();
261 
262  //**************************************************
263  //Place the single crystal in the 2x2 volume 4 times
264  //**************************************************
265 
266  G4int j_idx, k_idx;
267  G4double x_cent, y_cent, z_cent, rot_x, rot_y, rot_z;
268  G4int MappingIndex;
269 
270  j = 0;
271  while (j_cry > j)
272  {
273  MappingIndex = TwoByTwo[j][8];
274  if (MappingIndex == 1)
275  {
276  j_idx = TwoByTwo[j][0];
277  k_idx = TwoByTwo[j][1];
278  x_cent = TwoByTwo[j][2];
279  y_cent = TwoByTwo[j][3];
280  z_cent = TwoByTwo[j][4];
281  //rot_x = TwoByTwo[j][5];
282  //rot_y = TwoByTwo[j][6];
283  rot_z = TwoByTwo[j][7];
284 
285  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
286 
287  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
288  Rot->rotateX(0 * rad);
289  Rot->rotateY(0 * rad);
290  Rot->rotateZ(rot_z * rad);
291 
292  string crystal_name = _crystallogicnameprefix + "_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
293  int copyno = (j_idx << 16) + k_idx;
294  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot, Crystal_Center,
295  crystal_logic_small,
296  crystal_name,
297  Two_by_Two_logic,
298  0, copyno, OverlapCheck());
299  m_ActiveVolumeSet.insert(physvol);
300  j_idx = k_idx = 0;
301  x_cent = y_cent = z_cent = rot_x = rot_y = rot_z = 0.0;
302  }
303  j++;
304  }
305 
306  //*************************************************************
307  //Place the 2x2 volume in the 4x4 volume 4 times, with rotation
308  //*************************************************************
309 
310  j = 0;
311  while (j_cry > j)
312  {
313  MappingIndex = TwoByTwo[j][8];
314  if (MappingIndex == 4)
315  {
316  j_idx = TwoByTwo[j][0];
317  k_idx = TwoByTwo[j][1];
318  x_cent = TwoByTwo[j][2];
319  y_cent = TwoByTwo[j][3];
320  z_cent = TwoByTwo[j][4];
321  rot_x = TwoByTwo[j][5];
322  rot_y = TwoByTwo[j][6];
323  rot_z = TwoByTwo[j][7];
324 
325  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
326 
327  G4RotationMatrix *Rot = new G4RotationMatrix();
328  Rot->rotateX(rot_x * rad);
329  Rot->rotateY(rot_y * rad);
330  Rot->rotateZ(0 * rad);
331 
332  string Two_by_Two_name = "TwoByTwo_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
333 
334  int copyno = (j_idx << 16) + k_idx;
335  new G4PVPlacement(Rot, Crystal_Center,
336  Two_by_Two_logic,
337  Two_by_Two_name,
338  crystal_logic,
339  0, copyno, OverlapCheck());
340 
341  j_idx = k_idx = 0;
342  x_cent = y_cent = z_cent = rot_x = rot_y = 0.0;
343  }
344 
345  j++;
346  }
347 
348  //*****************************
349  //Create the carbon fiber shell
350  //*****************************
351 
352  //Create a hunk of carbon fiber the same size as the mother volume
353 
354  G4double dx1, dy1, dx2, dy2, dz_whole, carbon_fiber_adjust_width, carbon_fiber_adjust_length;
355 
356  GetCarbonFiberAdjustments(carbon_fiber_adjust_width, carbon_fiber_adjust_length);
357 
358  dx1 = _dx_front + carbon_fiber_adjust_width;
359  dy1 = dx1;
360  dx2 = _dx_back - carbon_fiber_adjust_width;
361  dy2 = dx2;
362  dz_whole = _dz_crystal - carbon_fiber_adjust_length; //Carbon fiber shell should be slightly shorter than whole volume
363 
364  G4VSolid *Carbon_hunk_solid = new G4Trd(G4String("Carbon_hunk_solid"),
365  dx1, //Half length on the small face in x
366  dx2, //Half length on the large face in x
367  dy1, //Half length on the small face in y
368  dy2, //Half length on the large face in y
369  dz_whole); //Half length in z
370 
371  //Use subtraction solid to remove the crystal volumes from the carbon fiber hunk
372 
373  //----------------------------------------------------------------------------------------------------
374  //First 2x2 crystal
375 
376  j = 0;
377  G4int counter = 0;
378  while (j_cry > counter)
379  {
380  MappingIndex = TwoByTwo[j][8];
381  if (MappingIndex != 4) j++;
382  counter++;
383  }
384 
385  x_cent = TwoByTwo[j][2];
386  y_cent = TwoByTwo[j][3];
387  z_cent = TwoByTwo[j][4];
388  rot_x = TwoByTwo[j][5];
389  rot_y = TwoByTwo[j][6];
390  //rot_z = TwoByTwo[j][7];
391 
392  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
393 
394  G4RotationMatrix *Rot_1 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
395  Rot_1->rotateX(rot_x * rad);
396  Rot_1->rotateY(rot_y * rad);
397  Rot_1->rotateZ(0 * rad);
398 
399  G4SubtractionSolid *Carbon_Shell_1 = new G4SubtractionSolid(G4String("Carbon_Shell_1"),
400  Carbon_hunk_solid,
401  Two_by_Two_solid,
402  Rot_1,
403  Crystal_Center);
404  j++;
405 
406  //----------------------------------------------------------------------------------------------------
407  //Second 2x2 crystal
408 
409  x_cent = TwoByTwo[j][2];
410  y_cent = TwoByTwo[j][3];
411  z_cent = TwoByTwo[j][4];
412  rot_x = TwoByTwo[j][5];
413  rot_y = TwoByTwo[j][6];
414  //rot_z = TwoByTwo[j][7];
415 
416  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
417 
418  G4RotationMatrix *Rot_2 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
419  Rot_2->rotateX(rot_x * rad);
420  Rot_2->rotateY(rot_y * rad);
421  Rot_2->rotateZ(0 * rad);
422 
423  G4SubtractionSolid *Carbon_Shell_2 = new G4SubtractionSolid(G4String("Carbon_Shell_2"),
424  Carbon_Shell_1,
425  Two_by_Two_solid,
426  Rot_2,
427  Crystal_Center);
428  j++;
429 
430  //----------------------------------------------------------------------------------------------------
431  //Third 2x2 crystal
432 
433  x_cent = TwoByTwo[j][2];
434  y_cent = TwoByTwo[j][3];
435  z_cent = TwoByTwo[j][4];
436  rot_x = TwoByTwo[j][5];
437  rot_y = TwoByTwo[j][6];
438  //rot_z = TwoByTwo[j][7];
439 
440  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
441 
442  G4RotationMatrix *Rot_3 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
443  Rot_3->rotateX(rot_x * rad);
444  Rot_3->rotateY(rot_y * rad);
445  Rot_3->rotateZ(0 * rad);
446 
447  G4SubtractionSolid *Carbon_Shell_3 = new G4SubtractionSolid(G4String("Carbon_Shell_3"),
448  Carbon_Shell_2,
449  Two_by_Two_solid,
450  Rot_3,
451  Crystal_Center);
452  j++;
453 
454  //----------------------------------------------------------------------------------------------------
455  //Final 2x2 crystal
456 
457  x_cent = TwoByTwo[j][2];
458  y_cent = TwoByTwo[j][3];
459  z_cent = TwoByTwo[j][4];
460  rot_x = TwoByTwo[j][5];
461  rot_y = TwoByTwo[j][6];
462  //rot_z = TwoByTwo[j][7];
463 
464  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
465 
466  G4RotationMatrix *Rot_4 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
467  Rot_4->rotateX(rot_x * rad);
468  Rot_4->rotateY(rot_y * rad);
469  Rot_4->rotateZ(0 * rad);
470 
471  G4SubtractionSolid *Carbon_Shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
472  Carbon_Shell_3,
473  Two_by_Two_solid,
474  Rot_4,
475  Crystal_Center);
476  j++;
477 
478  //----------------------------------------------------------------------------------------------------
479 
480  //Create logical volume with the subtracted solid, made from carbon fiber material defined earlier
481 
482  G4LogicalVolume *Carbon_Shell_logic = new G4LogicalVolume(Carbon_Shell_solid,
483  GetCarbonFiber(),
484  "Carbon_Fiber_logic",
485  0, 0, 0);
486 
487  GetDisplayAction()->AddVolume(Carbon_Shell_logic, "CarbonShell");
488 
489  //*************************************
490  //Place the carbon fiber shell at (0,0)
491  //*************************************
492 
493  G4ThreeVector Carbon_Center = G4ThreeVector(0 * mm, 0 * mm, 0 * mm);
494 
495  G4RotationMatrix *Rot_5 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
496  Rot_5->rotateX(0 * rad);
497  Rot_5->rotateY(0 * rad);
498  Rot_5->rotateZ(0 * rad);
499 
500  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot_5, Carbon_Center,
501  Carbon_Shell_logic,
502  "Carbon_Fiber_Shell",
503  crystal_logic,
504  0, 0, OverlapCheck());
505 
506  m_PassiveVolumeSet.insert(physvol);
507  //***********************************
508  //All done! Return to parent function
509  //***********************************
510 
511  return 0;
512 }
513 
514 int PHG4ProjCrystalCalorimeterDetector::FillSpecialUnit(G4LogicalVolume *crystal_logic, G4int ident)
515 {
516  //*************************************
517  //**********Define Materials***********
518  //*************************************
519 
520  //Crystal Material (Default is Lead Tungstate)
521  G4Material *material_crystal = GetDetectorMaterial(GetParams()->get_string_param("material"));
522 
523  recoConsts *rc = recoConsts::instance();
524  G4Material *WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
525 
526  //*************************************
527  //**********Build First Crystal********
528  //*************************************
529 
530  //Crystal Dimensions determined by the _dx_front, with various gaps and carbon fiber widths subtracted out
531 
532  //Parameters of the spacing given by PANDA document arXiv:0810.1216v1 Fig. 7.25
533  G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
534  GetCarbonFiberSpacing(carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals);
535 
536  //Crystal Dimensions
537  G4double dx_front_small = (_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the front crystal face
538  G4double dy_front_small = (_dy_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full height of the front crystal face
539  G4double dx_back_small = (_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the back crystal face
540  G4double dy_back_small = (_dy_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full height of the back crystal face
541  G4double dz = _dz_crystal;
542 
543  //Vertices of the primary, irregularly shaped crystal put into a vector
544  std::vector<G4TwoVector> vertices;
545  vertices.push_back(G4TwoVector(0, 0));
546  vertices.push_back(G4TwoVector(0, dy_front_small));
547  vertices.push_back(G4TwoVector(dx_front_small, dy_front_small));
548  vertices.push_back(G4TwoVector(dx_front_small, 0));
549  vertices.push_back(G4TwoVector(0, 0));
550  vertices.push_back(G4TwoVector(0, dy_back_small));
551  vertices.push_back(G4TwoVector(dx_back_small, dy_back_small));
552  vertices.push_back(G4TwoVector(dx_back_small, 0));
553 
554  //Create the primary, irregularly shaped crystal
555  G4VSolid *crystal_solid_small = new G4GenericTrap(G4String("eEcal_crystal"),
556  dz,
557  vertices);
558 
559  G4LogicalVolume *crystal_logic_small = new G4LogicalVolume(crystal_solid_small,
560  material_crystal,
561  "eEcal_crystal",
562  0, 0, 0);
563 
564  GetDisplayAction()->AddVolume(crystal_logic_small, "Crystal");
565 
566  //****************************************************
567  //Build the solid/logical volume for the 2 x 2 crystal
568  //****************************************************
569 
570  G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
571  G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
572  G4double TwoByTwo_dy1 = TwoByTwo_dx1;
573  G4double TwoByTwo_dy2 = TwoByTwo_dx2;
574  G4double TwoByTwo_dz = _dz_crystal;
575 
576  G4VSolid *Two_by_Two_solid = new G4Trd(G4String("Two_by_Two_solid"),
577  TwoByTwo_dx1, //Half length on the small face in x
578  TwoByTwo_dx2, //Half length on the large face in x
579  TwoByTwo_dy1, //Half length on the small face in y
580  TwoByTwo_dy2, //Half length on the large face in y
581  TwoByTwo_dz); //Half length in z
582 
583  G4LogicalVolume *Two_by_Two_logic = new G4LogicalVolume(Two_by_Two_solid,
584  WorldMaterial,
585  "2_by_2_unit",
586  0, 0, 0);
587 
588  GetDisplayAction()->AddVolume(Two_by_Two_logic, "TwoByTwo");
589 
590  //*************************************************************
591  //Read in mapping file for a single 2 x 2 block and 4 x 4 block
592  //*************************************************************
593 
594  //The first four lines of the data file refer to the 2x2 block, and the last four lines refer to the mapping of the 4x4 block
595 
596  const int NumberOfIndices = 9; //Number of indices in mapping file for 4x4 block
597 
598  //Find the number of lines in the file, make and fill a NumberOfLines by NumberOfIndices matrix with contents of data file
599  int NumberOfLines = 0;
600  ifstream in(GetParams()->get_string_param("mapping4x4"));
601  if (!in.is_open())
602  {
603  cout << endl
604  << "*******************************************************************" << endl;
605  cout << "ERROR in 2 by 2 crystal mapping ";
606  cout << "Failed to open " << GetParams()->get_string_param("mapping4x4") << " --- Exiting program." << endl;
607  cout << "*******************************************************************" << endl
608  << endl;
609  gSystem->Exit(1);
610  }
611  std::string unused;
612  while (std::getline(in, unused))
613  {
614  ++NumberOfLines;
615  }
616  in.close();
617 
618  in.open(GetParams()->get_string_param("mapping4x4"));
619  G4int j_cry = NumberOfLines;
620  G4int k_cry = NumberOfIndices;
621 
622  double TwoByTwo[j_cry][k_cry];
623 
624  G4int j = 0;
625  G4int k = 0;
626 
627  while (j_cry > j)
628  {
629  while (k_cry > k)
630  {
631  in >> TwoByTwo[j][k];
632  k++;
633  }
634  j++;
635  k = 0;
636  }
637  in.close();
638  //**************************************************
639  //Place the single crystal in the 2x2 volume 4 times
640  //**************************************************
641 
642  G4int j_idx, k_idx;
643  G4double x_cent, y_cent, z_cent, rot_x, rot_y, rot_z;
644  G4int MappingIndex;
645 
646  j = 0;
647  while (j_cry > j)
648  {
649  MappingIndex = TwoByTwo[j][8];
650  if (MappingIndex == 1)
651  {
652  j_idx = TwoByTwo[j][0];
653  k_idx = TwoByTwo[j][1];
654  x_cent = TwoByTwo[j][2];
655  y_cent = TwoByTwo[j][3];
656  z_cent = TwoByTwo[j][4];
657  rot_x = TwoByTwo[j][5];
658  rot_y = TwoByTwo[j][6];
659  rot_z = TwoByTwo[j][7];
660 
661  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
662 
663  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
664  Rot->rotateX(0 * rad);
665  Rot->rotateY(0 * rad);
666  Rot->rotateZ(rot_z * rad);
667 
668  string crystal_name = _crystallogicnameprefix + "_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
669  int copyno = (j_idx << 16) + k_idx;
670 
671  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot, Crystal_Center,
672  crystal_logic_small,
673  crystal_name,
674  Two_by_Two_logic,
675  0, copyno, OverlapCheck());
676  m_ActiveVolumeSet.insert(physvol);
677 
678  j_idx = k_idx = 0;
679  x_cent = y_cent = z_cent = rot_z = 0.0;
680  }
681  j++;
682  }
683 
684  //****************************
685  //Place the 2x2 Crystal Blocks
686  //****************************
687 
688  j = 0;
689  while (j_cry > j)
690  {
691  j_idx = TwoByTwo[j][0];
692  k_idx = TwoByTwo[j][1];
693  x_cent = TwoByTwo[j][2];
694  y_cent = TwoByTwo[j][3];
695  z_cent = TwoByTwo[j][4];
696  rot_x = TwoByTwo[j][5];
697  rot_y = TwoByTwo[j][6];
698  rot_z = TwoByTwo[j][7];
699  MappingIndex = TwoByTwo[j][8];
700 
701  if (ident == 12)
702  {
703  G4ThreeVector Crystal_Center = G4ThreeVector(0 * mm, 0 * mm, 0 * mm);
704 
705  G4RotationMatrix *Rot = new G4RotationMatrix();
706  Rot->rotateX(0 * rad);
707  Rot->rotateY(0 * rad);
708  Rot->rotateZ(0 * rad);
709 
710  string Two_by_Two_name = "TwoByTwo_j_0_k_0";
711  int copyno = 0;
712  new G4PVPlacement(Rot, Crystal_Center,
713  Two_by_Two_logic,
714  Two_by_Two_name,
715  crystal_logic,
716  0, copyno, OverlapCheck());
717  }
718  else if (ident == 22)
719  {
720  if (MappingIndex == 22)
721  {
722  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
723 
724  G4RotationMatrix *Rot = new G4RotationMatrix();
725  Rot->rotateX(rot_x * rad);
726  Rot->rotateY(0 * rad);
727  Rot->rotateZ(0 * rad);
728 
729  string Two_by_Two_name = "TwoByTwo_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
730  int copyno = (j_idx << 16) + k_idx;
731  new G4PVPlacement(Rot, Crystal_Center,
732  Two_by_Two_logic,
733  Two_by_Two_name,
734  crystal_logic,
735  0, copyno, OverlapCheck());
736  }
737  }
738  else if (ident == 32)
739  {
740  if (MappingIndex == 32)
741  {
742  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
743 
744  G4RotationMatrix *Rot = new G4RotationMatrix();
745  Rot->rotateX(rot_x * rad);
746  Rot->rotateY(rot_y * rad);
747  Rot->rotateZ(0 * rad);
748 
749  string Two_by_Two_name = "TwoByTwo_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
750  int copyno = (j_idx << 16) + k_idx;
751  new G4PVPlacement(Rot, Crystal_Center,
752  Two_by_Two_logic,
753  Two_by_Two_name,
754  crystal_logic,
755  0, copyno, OverlapCheck());
756  }
757  }
758  else
759  {
760  cerr << endl
761  << "Invalid Mapping Type: " << ident << endl;
762  return -1;
763  }
764 
765  j_idx = k_idx = 0;
766  x_cent = y_cent = z_cent = rot_x = rot_y = rot_z = 0.0;
767  j++;
768  }
769 
770  //*******************
771  //Create Carbon Fiber
772  //*******************
773 
774  G4double dx1, dx2, dy1, dy2;
775  GetCrystalSize(dx1, dy1, dx2, dy2, dz);
776 
777  if (ident == 12)
778  {
779  G4VSolid *Carbon_hunk_solid = new G4Trd(G4String("Carbon_hunk_solid"),
780  (dx1 / 2.00 + 0.0209292929), //Half length on the small face in x
781  (dx2 / 2.00 - 0.0209292929), //Half length on the large face in x
782  (dy1 / 2.00 + 0.0209292929), //Half length on the small face in y
783  (dy2 / 2.0 - 0.0209292929), //Half length on the large face in y
784  (dz - 2 * mm));
785 
786  x_cent = 0 * mm;
787  y_cent = 0 * mm;
788  z_cent = 0 * mm;
789  rot_x = 0 * mm;
790  rot_y = 0 * mm;
791  rot_z = 0 * mm;
792 
793  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
794 
795  G4RotationMatrix *Rot_1 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
796  Rot_1->rotateX(rot_x * rad);
797  Rot_1->rotateY(rot_y * rad);
798  Rot_1->rotateZ(rot_z * rad);
799 
800  G4SubtractionSolid *Carbon_Shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
801  Carbon_hunk_solid,
802  Two_by_Two_solid,
803  Rot_1,
804  Crystal_Center);
805 
806  G4LogicalVolume *Carbon_Shell_logic = new G4LogicalVolume(Carbon_Shell_solid,
807  GetCarbonFiber(),
808  "Carbon_Fiber_logic",
809  0, 0, 0);
810 
811  GetDisplayAction()->AddVolume(Carbon_Shell_logic, "CarbonShell");
812 
813  //*************************************
814  //Place the carbon fiber shell at (0,0)
815  //*************************************
816 
817  G4ThreeVector Carbon_Center = G4ThreeVector(0 * mm, 0 * mm, 0 * mm);
818 
819  G4RotationMatrix *Rot_5 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
820  Rot_5->rotateX(0 * rad);
821  Rot_5->rotateY(0 * rad);
822  Rot_5->rotateZ(0 * rad);
823 
824  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot_5, Carbon_Center,
825  Carbon_Shell_logic,
826  "Carbon_Fiber_Shell",
827  crystal_logic,
828  0, 0, OverlapCheck());
829  m_PassiveVolumeSet.insert(physvol);
830  }
831  else if (ident == 22)
832  {
833  G4double carbon_fiber_adjust_width; //Because the crystals are slightly angled, the carbon fiber needs to be shortened
834  G4double carbon_fiber_adjust_length; // from the mother volume (to prevent clipping) by this amount.
835  GetCarbonFiberAdjustments(carbon_fiber_adjust_width, carbon_fiber_adjust_length);
836  G4double x_adjust = 0.0519558696 * mm;
837 
838  G4VSolid *Carbon_hunk_solid = new G4Trd(G4String("Carbon_hunk_solid"),
839  (dx1 / 2.0 + (x_adjust)), //Half length on the small face in x
840  (dx2 / 2.0 - (x_adjust)), //Half length on the large face in x
841  (dy1 + carbon_fiber_adjust_width), //Half length on the small face in y
842  (dy2 - carbon_fiber_adjust_width), //Half length on the large face in y
843  (dz - carbon_fiber_adjust_length)); //Half length in z
844 
845  j = 0;
846  G4int counter = 0;
847  while (j_cry > counter)
848  {
849  MappingIndex = TwoByTwo[j][8];
850  if (MappingIndex != 22) j++;
851  counter++;
852  }
853 
854  //First Hole
855 
856  x_cent = TwoByTwo[j][2];
857  y_cent = TwoByTwo[j][3];
858  z_cent = TwoByTwo[j][4];
859  rot_x = TwoByTwo[j][5];
860  rot_y = TwoByTwo[j][6];
861  rot_z = TwoByTwo[j][7];
862 
863  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
864 
865  G4RotationMatrix *Rot_1 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
866  Rot_1->rotateX(rot_x * rad);
867  Rot_1->rotateY(0 * rad);
868  Rot_1->rotateZ(0 * rad);
869 
870  G4SubtractionSolid *Carbon_Shell_1 = new G4SubtractionSolid(G4String("Carbon_Shell_1"),
871  Carbon_hunk_solid,
872  Two_by_Two_solid,
873  Rot_1,
874  Crystal_Center);
875 
876  j++;
877 
878  //Second Hole
879 
880  x_cent = TwoByTwo[j][2];
881  y_cent = TwoByTwo[j][3];
882  z_cent = TwoByTwo[j][4];
883  rot_x = TwoByTwo[j][5];
884  //rot_y = TwoByTwo[j][6];
885  //rot_z = TwoByTwo[j][7];
886 
887  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
888 
889  G4RotationMatrix *Rot_2 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
890  Rot_2->rotateX(rot_x * rad);
891  Rot_2->rotateY(0 * rad);
892  Rot_2->rotateZ(0 * rad);
893 
894  G4SubtractionSolid *Carbon_Shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
895  Carbon_Shell_1,
896  Two_by_Two_solid,
897  Rot_2,
898  Crystal_Center);
899 
900  G4LogicalVolume *Carbon_Shell_logic = new G4LogicalVolume(Carbon_Shell_solid,
901  GetCarbonFiber(),
902  "Carbon_Fiber_logic",
903  0, 0, 0);
904 
905  GetDisplayAction()->AddVolume(Carbon_Shell_logic, "CarbonShell");
906 
907  //*************************************
908  //Place the carbon fiber shell at (0,0)
909  //*************************************
910 
911  G4ThreeVector Carbon_Center = G4ThreeVector(0 * mm, 0 * mm, 0 * mm);
912 
913  G4RotationMatrix *Rot_5 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
914  Rot_5->rotateX(0 * rad);
915  Rot_5->rotateY(0 * rad);
916  Rot_5->rotateZ(0 * rad);
917 
918  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot_5, Carbon_Center,
919  Carbon_Shell_logic,
920  "Carbon_Fiber_Shell",
921  crystal_logic,
922  0, 0, OverlapCheck());
923  m_PassiveVolumeSet.insert(physvol);
924  }
925  else if (ident == 32)
926  {
927  x_cent = 22.6036363636 + 0.18000 * 2.0;
928  y_cent = x_cent;
929  z_cent = 0.00;
930  G4double rot_x = 0.020928529;
931  G4double rot_y = -1.0 * rot_x;
932 
933  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
934 
935  G4RotationMatrix *Rot_1 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
936  Rot_1->rotateX(rot_x * rad);
937  Rot_1->rotateY(rot_y * rad);
938  Rot_1->rotateZ(0 * rad);
939 
940  G4VSolid *FourByFour_hunk_solid = new G4Trd(G4String("4x4_hunk_solid"),
941  dx1, //Half length on the small face in x
942  dx2, //Half length on the large face in x
943  dy1, //Half length on the small face in y
944  dy2, //Half length on the large face in y
945  dz); //Half length in z
946 
947  //Parameters of the spacing given by PANDA document arXiv:0810.1216v1 Fig. 7.25
948  G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
949  GetCarbonFiberSpacing(carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals);
950 
951  //Crystal Dimensions
952  G4double dx_front_small = (_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the front crystal face
953  //G4double dy_front_small = ( _dy_front - (2.0 * carbon_fiber_width ) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals ) / 2.0; //Full height of the front crystal face
954  G4double dx_back_small = (_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the back crystal face
955  //G4double dy_back_small = (_dy_back - (2.0 * carbon_fiber_width ) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals ) / 2.0; //Full height of the back crystal face
956 
957  G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
958  G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
959  G4double TwoByTwo_dy1 = TwoByTwo_dx1;
960  G4double TwoByTwo_dy2 = TwoByTwo_dx2;
961  G4double TwoByTwo_dz = _dz_crystal;
962 
963  G4VSolid *Two_by_Two_solid = new G4Trd(G4String("Two_by_Two_solid"),
964  TwoByTwo_dx1, //Half length on the small face in x
965  TwoByTwo_dx2, //Half length on the large face in x
966  TwoByTwo_dy1, //Half length on the small face in y
967  TwoByTwo_dy2, //Half length on the large face in y
968  TwoByTwo_dz + 2 * mm); //Half length in z
969 
970  G4SubtractionSolid *Carbon_hunk_solid = new G4SubtractionSolid(G4String("Carbon_hunk_solid"),
971  FourByFour_hunk_solid,
972  Two_by_Two_solid,
973  Rot_1,
974  Crystal_Center);
975 
976  //----------------------------------------------------------------------------------------------------
977  //First 2x2 crystal hole
978 
979  j = 0;
980  G4int counter = 0;
981  while (j_cry > counter)
982  {
983  MappingIndex = TwoByTwo[j][8];
984  if (MappingIndex != 32) j++;
985  counter++;
986  }
987 
988  x_cent = TwoByTwo[j][2];
989  y_cent = TwoByTwo[j][3];
990  z_cent = TwoByTwo[j][4];
991  rot_x = TwoByTwo[j][5];
992  rot_y = TwoByTwo[j][6];
993  rot_z = TwoByTwo[j][7];
994 
995  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
996 
997  G4RotationMatrix *Rot_6 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
998  Rot_6->rotateX(rot_x * rad);
999  Rot_6->rotateY(rot_y * rad);
1000  Rot_6->rotateZ(0 * rad);
1001 
1002  G4SubtractionSolid *Carbon_Shell_1 = new G4SubtractionSolid(G4String("Carbon_Shell_1"),
1003  Carbon_hunk_solid,
1004  Two_by_Two_solid,
1005  Rot_6,
1006  Crystal_Center);
1007  j++;
1008 
1009  //----------------------------------------------------------------------------------------------------
1010  //Second 2x2 crystal hole
1011 
1012  x_cent = TwoByTwo[j][2];
1013  y_cent = TwoByTwo[j][3];
1014  z_cent = TwoByTwo[j][4];
1015  rot_x = TwoByTwo[j][5];
1016  rot_y = TwoByTwo[j][6];
1017  //rot_z = TwoByTwo[j][7];
1018 
1019  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1020 
1021  G4RotationMatrix *Rot_2 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1022  Rot_2->rotateX(rot_x * rad);
1023  Rot_2->rotateY(rot_y * rad);
1024  Rot_2->rotateZ(0 * rad);
1025 
1026  G4SubtractionSolid *Carbon_Shell_2 = new G4SubtractionSolid(G4String("Carbon_Shell_2"),
1027  Carbon_Shell_1,
1028  Two_by_Two_solid,
1029  Rot_2,
1030  Crystal_Center);
1031  j++;
1032 
1033  //----------------------------------------------------------------------------------------------------
1034  //Third 2x2 crystal hole
1035 
1036  x_cent = TwoByTwo[j][2];
1037  y_cent = TwoByTwo[j][3];
1038  z_cent = TwoByTwo[j][4];
1039  rot_x = TwoByTwo[j][5];
1040  rot_y = TwoByTwo[j][6];
1041  //rot_z = TwoByTwo[j][7];
1042 
1043  Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1044 
1045  G4RotationMatrix *Rot_3 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1046  Rot_3->rotateX(rot_x * rad);
1047  Rot_3->rotateY(rot_y * rad);
1048  Rot_3->rotateZ(0 * rad);
1049 
1050  G4SubtractionSolid *Carbon_Shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
1051  Carbon_Shell_2,
1052  Two_by_Two_solid,
1053  Rot_3,
1054  Crystal_Center);
1055 
1056  G4LogicalVolume *Carbon_Shell_logic = new G4LogicalVolume(Carbon_Shell_solid,
1057  GetCarbonFiber(),
1058  "Carbon_Fiber_logic",
1059  0, 0, 0);
1060 
1061  GetDisplayAction()->AddVolume(Carbon_Shell_logic, "CarbonShell");
1062 
1063  //*************************************
1064  //Place the carbon fiber shell at (0,0)
1065  //*************************************
1066 
1067  G4ThreeVector Carbon_Center = G4ThreeVector(0 * mm, 0 * mm, 0 * mm);
1068 
1069  G4RotationMatrix *Rot_5 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1070  Rot_5->rotateX(0 * rad);
1071  Rot_5->rotateY(0 * rad);
1072  Rot_5->rotateZ(0 * rad);
1073 
1074  G4VPhysicalVolume *physvol = new G4PVPlacement(Rot_5, Carbon_Center,
1075  Carbon_Shell_logic,
1076  "Carbon_Fiber_Shell",
1077  crystal_logic,
1078  0, 0, OverlapCheck());
1079  m_PassiveVolumeSet.insert(physvol);
1080  }
1081  else
1082  {
1083  cerr << endl
1084  << "This is an error message which should never be shown. You may have disabled an important portion of this code upsream. Sorry :)" << endl;
1085  return -1;
1086  }
1087 
1088  return 0;
1089 }
1090 
1091 //_______________________________________________________________________
1092 int PHG4ProjCrystalCalorimeterDetector::ConstructProjectiveCrystals(G4LogicalVolume *ecalenvelope)
1093 {
1094  G4int NumberOfLines; //Number of crystals to be created.
1095  const G4int NumberOfIndices = 9; //Different dimensions needed for crystal placement
1096  const string FileName = GetParams()->get_string_param("mappingtower"); //File in which crystal positions are stored
1097 
1098  G4int j_cry, k_cry; //Indices for matrix
1099  G4int j_idx, k_idx; //Indices of each crstals
1100  G4double x_cent, y_cent, z_cent, r_theta, r_phi; //Coordinates of crystal in [x,y,z,theta,phi]
1101 
1102  G4double dx1; //Half of the extent of the front face of the trapezoid in x
1103  G4double dx2; //Half of the extent of the back face of the trapezoid in x
1104  G4double dy1; //Half of the extent of the front face of the trapezoid in y
1105  G4double dy2; //Half of the extent of the back face of the trapezoid in y
1106  G4double dz; //Half of the extent of the crystal in z
1107 
1108  GetCrystalSize(dx1, dy1, dx2, dy2, dz); //Fill crystal dimensions with function PHG4ProjCrystalCalorimeterDetector::CrystalDimensions
1109 
1110  //Create single crystal
1111 
1112  G4VSolid *crystal_solid = new G4Trd(G4String("eEcal_crystal"),
1113  dx1, //Half length on the small face in x
1114  dx2, //Half length on the large face in x
1115  dy1, //Half length on the small face in y
1116  dy2, //Half length on the large face in y
1117  dz); //Half length in z
1118 
1119  recoConsts *rc = recoConsts::instance();
1120  G4Material *WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
1121 
1122  G4LogicalVolume *crystal_logic = new G4LogicalVolume(crystal_solid,
1123  WorldMaterial,
1124  "eEcal_crystal_unit",
1125  0, 0, 0);
1126 
1127  //Create the single 2x2 unit (MappingIndex = 12)
1128 
1129  G4VSolid *twelve_solid = new G4Trd(G4String("12_solid"),
1130  (dx1 / 2.0), //Half length on the small face in x
1131  (dx2 / 2.0), //Half length on the large face in x
1132  (dy1 / 2.0), //Half length on the small face in y
1133  (dy2 / 2.0), //Half length on the large face in y
1134  dz); //Half length in z
1135 
1136  G4LogicalVolume *twelve_logic = new G4LogicalVolume(twelve_solid,
1137  WorldMaterial,
1138  "12_unit",
1139  0, 0, 0);
1140 
1141  //Create the double 2x2 unit (i-shaped, MappingIndex = 22)
1142 
1143  G4VSolid *twentytwo_solid = new G4Trd(G4String("22_solid"),
1144  (dx1 / 2.0), //Half length on the small face in x
1145  (dx2 / 2.0), //Half length on the large face in x
1146  dy1, //Half length on the small face in y
1147  dy2, //Half length on the large face in y
1148  dz); //Half length in z
1149 
1150  G4LogicalVolume *twentytwo_logic = new G4LogicalVolume(twentytwo_solid,
1151  WorldMaterial,
1152  "22_unit",
1153  0, 0, 0);
1154 
1155  //Create the triple 2x2 unit (L-shaped, MappingIndex = 32)
1156 
1157  x_cent = 22.6036363636 + 0.18000 * 2.0;
1158  y_cent = x_cent;
1159  z_cent = 0.00;
1160  G4double rot_x = 0.020928529;
1161  G4double rot_y = -1.0 * rot_x;
1162  G4double rot_z = 0.00;
1163 
1164  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1165 
1166  G4RotationMatrix *Rot_1 = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1167  Rot_1->rotateX(rot_x * rad);
1168  Rot_1->rotateY(rot_y * rad);
1169  Rot_1->rotateZ(0 * rad);
1170 
1171  G4VSolid *FourByFour_hunk_solid = new G4Trd(G4String("4x4_hunk_solid"),
1172  dx1, //Half length on the small face in x
1173  dx2, //Half length on the large face in x
1174  dy1, //Half length on the small face in y
1175  dy2, //Half length on the large face in y
1176  dz); //Half length in z
1177 
1178  //Parameters of the spacing given by PANDA document arXiv:0810.1216v1 Fig. 7.25
1179  G4double carbon_fiber_width = 0.18 * mm; //Width of the carbon fiber which surrounds the crystal
1180  G4double air_gap_carbon_fiber = 0.24 * mm; //Air gap between crystal and the carbon fiber
1181  G4double air_gap_crystals = 0.60 * mm; //Air gap between crystal and crystal
1182 
1183  //Crystal Dimensions
1184  G4double dx_front_small = (_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the front crystal face
1185  // G4double dy_front_small = ( _dy_front - (2.0 * carbon_fiber_width ) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals ) / 2.0; //Full height of the front crystal face
1186  G4double dx_back_small = (_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0; //Full width of the back crystal face
1187  // G4double dy_back_small = (_dy_back - (2.0 * carbon_fiber_width ) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals ) / 2.0; //Full height of the back crystal face
1188 
1189  G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
1190  G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
1191  G4double TwoByTwo_dy1 = TwoByTwo_dx1;
1192  G4double TwoByTwo_dy2 = TwoByTwo_dx2;
1193  G4double TwoByTwo_dz = _dz_crystal;
1194 
1195  G4VSolid *Two_by_Two_solid = new G4Trd(G4String("Two_by_Two_solid"),
1196  TwoByTwo_dx1, //Half length on the small face in x
1197  TwoByTwo_dx2, //Half length on the large face in x
1198  TwoByTwo_dy1, //Half length on the small face in y
1199  TwoByTwo_dy2, //Half length on the large face in y
1200  TwoByTwo_dz + 5 * mm); //Half length in z
1201 
1202  G4SubtractionSolid *thirtytwo_solid = new G4SubtractionSolid(G4String("32_solid"),
1203  FourByFour_hunk_solid,
1204  Two_by_Two_solid,
1205  Rot_1,
1206  Crystal_Center);
1207 
1208  G4LogicalVolume *thirtytwo_logic = new G4LogicalVolume(thirtytwo_solid,
1209  WorldMaterial,
1210  "32_unit",
1211  0, 0, 0);
1212 
1213  GetDisplayAction()->AddVolume(crystal_logic, "Invisible");
1214  GetDisplayAction()->AddVolume(twelve_logic, "Invisible");
1215  GetDisplayAction()->AddVolume(twentytwo_logic, "Invisible");
1216  GetDisplayAction()->AddVolume(thirtytwo_logic, "Invisible");
1217 
1218  Fill4x4Unit(crystal_logic);
1219  FillSpecialUnit(twelve_logic, 12);
1220  FillSpecialUnit(twentytwo_logic, 22);
1221  FillSpecialUnit(thirtytwo_logic, 32);
1222 
1223  ostringstream name;
1224 
1225  ifstream datafile;
1226 
1227  if (!datafile.is_open())
1228  {
1229  datafile.open(GetParams()->get_string_param("mappingtower"));
1230  if (!datafile)
1231  {
1232  cout << endl
1233  << "*******************************************************************" << endl;
1234  cout << "ERROR: Failed to open " << GetParams()->get_string_param("mappingtower") << " --- Exiting program." << endl;
1235  cout << "*******************************************************************" << endl
1236  << endl;
1237  gSystem->Exit(1);
1238  }
1239  }
1240 
1241  //Determine the number of crystals to be created
1242  NumberOfLines = 0;
1243  ifstream in(GetParams()->get_string_param("mappingtower"));
1244  std::string unused;
1245  while (std::getline(in, unused))
1246  ++NumberOfLines;
1247 
1248  j_cry = NumberOfLines; // = Number of Crystals
1249  k_cry = NumberOfIndices; // = j, k, x, y, z, alpha, beta.
1250 
1251  double Crystals[j_cry][k_cry];
1252 
1253  G4int j = 0;
1254  G4int k = 0;
1255 
1256  //Fill matrix with the data from the mapping file
1257  while (j_cry > j)
1258  {
1259  while (k_cry > k)
1260  {
1261  datafile >> Crystals[j][k];
1262  k++;
1263  }
1264  j++;
1265  k = 0;
1266  }
1267 
1268  // Find the maximum k index
1269  G4int k_max = 0;
1270  G4int MappingIndex;
1271  j = 0;
1272  while (j_cry > j)
1273  {
1274  MappingIndex = Crystals[j][8];
1275  if (Crystals[j][1] > k_max) k_max = Crystals[j][1];
1276  j++;
1277  }
1278 
1279  //Second Quadrant
1280  j = 0;
1281  while (j_cry > j)
1282  {
1283  MappingIndex = Crystals[j][8];
1284 
1285  j_idx = Crystals[j][0];
1286  k_idx = Crystals[j][1];
1287  x_cent = Crystals[j][2] - GetParams()->get_double_param("place_x") * cm;
1288  y_cent = Crystals[j][3] - GetParams()->get_double_param("place_y") * cm;
1289  z_cent = Crystals[j][4] - GetParams()->get_double_param("place_z") * cm; //Coordinate system refers to mother volume, have to subtract out its position in the actual xyz-space
1290  r_theta = Crystals[j][5]; //Rotation in Horizontal
1291  r_phi = Crystals[j][6]; //Rotation in Vertical
1292  rot_z = Crystals[j][7];
1293 
1294  if (MappingIndex == 16)
1295  {
1296  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1297 
1298  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1299  Rot->rotateX(r_phi * rad);
1300  Rot->rotateY(r_theta * rad);
1301  Rot->rotateZ(0 * rad);
1302 
1303  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1304  int copyno = (j_idx << 16) + k_idx;
1305  new G4PVPlacement(Rot, Crystal_Center,
1306  crystal_logic,
1307  FourByFourName,
1308  ecalenvelope,
1309  0, copyno, OverlapCheck());
1310  }
1311  else if (MappingIndex == 32)
1312  {
1313  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1314 
1315  if (j_idx == (k_max - 1) / 2 && k_idx == (k_max + 1) / 2) rot_z = rot_z / 2.0;
1316 
1317  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1318  Rot->rotateX(r_phi * rad);
1319  Rot->rotateY(r_theta * rad);
1320  Rot->rotateZ(rot_z * rad);
1321 
1322  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1323  int copyno = (j_idx << 16) + k_idx;
1324 
1325  new G4PVPlacement(Rot, Crystal_Center,
1326  thirtytwo_logic,
1327  FourByFourName,
1328  ecalenvelope,
1329  0, copyno, OverlapCheck());
1330  }
1331  else if (MappingIndex == 22)
1332  {
1333  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1334 
1335  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1336  Rot->rotateX(r_phi * rad);
1337  Rot->rotateY(r_theta * rad);
1338  Rot->rotateZ(rot_z * rad);
1339 
1340  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1341  int copyno = (j_idx << 16) + k_idx;
1342  new G4PVPlacement(Rot, Crystal_Center,
1343  twentytwo_logic,
1344  FourByFourName,
1345  ecalenvelope,
1346  0, copyno, OverlapCheck());
1347  }
1348  else if (MappingIndex == 12)
1349  {
1350  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1351 
1352  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1353  Rot->rotateX(r_phi * rad);
1354  Rot->rotateY(r_theta * rad);
1355  Rot->rotateZ(rot_z * rad);
1356 
1357  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1358  int copyno = (j_idx << 16) + k_idx;
1359 
1360  new G4PVPlacement(Rot, Crystal_Center,
1361  twelve_logic,
1362  FourByFourName,
1363  ecalenvelope,
1364  0, copyno, OverlapCheck());
1365  }
1366 
1367  j++;
1368  }
1369 
1370  //First Quadrant
1371  j = 0;
1372  while (j_cry > j)
1373  {
1374  MappingIndex = Crystals[j][8];
1375 
1376  j_idx = k_max - Crystals[j][0];
1377  k_idx = Crystals[j][1];
1378  x_cent = -1.0 * (Crystals[j][2] - GetParams()->get_double_param("place_x") * cm);
1379  y_cent = Crystals[j][3] - GetParams()->get_double_param("place_y") * cm;
1380  z_cent = Crystals[j][4] - GetParams()->get_double_param("place_z") * cm;
1381  r_theta = -1.0 * Crystals[j][5];
1382  r_phi = Crystals[j][6];
1383  rot_z = Crystals[j][7] / (-4.0);
1384 
1385  if (MappingIndex == 16)
1386  {
1387  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1388 
1389  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1390  Rot->rotateX(r_phi * rad);
1391  Rot->rotateY(r_theta * rad);
1392  Rot->rotateZ(0 * rad);
1393 
1394  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1395  int copyno = (j_idx << 16) + k_idx;
1396 
1397  new G4PVPlacement(Rot, Crystal_Center,
1398  crystal_logic,
1399  FourByFourName,
1400  ecalenvelope,
1401  0, copyno, OverlapCheck());
1402  }
1403  else if (MappingIndex == 32)
1404  {
1405  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1406 
1407  if (j_idx == (k_max + 1) / 2 && k_idx == (k_max + 1) / 2) rot_z = rot_z * 3.0;
1408 
1409  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1410  Rot->rotateX(r_phi * rad);
1411  Rot->rotateY(r_theta * rad);
1412  Rot->rotateZ(rot_z * rad);
1413 
1414  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1415  int copyno = (j_idx << 16) + k_idx;
1416  new G4PVPlacement(Rot, Crystal_Center,
1417  thirtytwo_logic,
1418  FourByFourName,
1419  ecalenvelope,
1420  0, copyno, OverlapCheck());
1421  }
1422  else if (MappingIndex == 22)
1423  {
1424  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1425 
1426  rot_z = Crystals[j][7] * (-1.0);
1427 
1428  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1429  Rot->rotateX(r_phi * rad);
1430  Rot->rotateY(r_theta * rad);
1431  Rot->rotateZ(rot_z * rad);
1432 
1433  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1434  int copyno = (j_idx << 16) + k_idx;
1435 
1436  new G4PVPlacement(Rot, Crystal_Center,
1437  twentytwo_logic,
1438  FourByFourName,
1439  ecalenvelope,
1440  0, copyno, OverlapCheck());
1441  }
1442  else if (MappingIndex == 12)
1443  {
1444  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1445 
1446  rot_z = rot_z * -4.0000;
1447 
1448  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1449  Rot->rotateX(r_phi * rad);
1450  Rot->rotateY(r_theta * rad);
1451  Rot->rotateZ(rot_z * rad);
1452 
1453  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1454  int copyno = (j_idx << 16) + k_idx;
1455 
1456  new G4PVPlacement(Rot, Crystal_Center,
1457  twelve_logic,
1458  FourByFourName,
1459  ecalenvelope,
1460  0, copyno, OverlapCheck());
1461  }
1462 
1463  j++;
1464  }
1465 
1466  //Fourth Quadrant
1467  j = 0;
1468  while (j_cry > j)
1469  {
1470  MappingIndex = Crystals[j][8];
1471  j_idx = k_max - Crystals[j][0];
1472  k_idx = k_max - Crystals[j][1];
1473  x_cent = -1.0 * (Crystals[j][2] - GetParams()->get_double_param("place_x") * cm);
1474  y_cent = -1.0 * (Crystals[j][3] - GetParams()->get_double_param("place_y") * cm);
1475  z_cent = Crystals[j][4] - GetParams()->get_double_param("place_z") * cm;
1476  r_theta = -1.0 * Crystals[j][5];
1477  r_phi = -1.0 * Crystals[j][6];
1478  rot_z = Crystals[j][7] / 2.0;
1479 
1480  if (MappingIndex == 16)
1481  {
1482  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1483 
1484  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1485  Rot->rotateX(r_phi * rad);
1486  Rot->rotateY(r_theta * rad);
1487  Rot->rotateZ(0 * rad);
1488 
1489  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1490  int copyno = (j_idx << 16) + k_idx;
1491 
1492  new G4PVPlacement(Rot, Crystal_Center,
1493  crystal_logic,
1494  FourByFourName,
1495  ecalenvelope,
1496  0, copyno, OverlapCheck());
1497  }
1498  else if (MappingIndex == 32)
1499  {
1500  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1501 
1502  if (j_idx == (k_max + 1) / 2 && k_idx == (k_max - 1) / 2) rot_z = rot_z * -2.0;
1503 
1504  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1505  Rot->rotateX(r_phi * rad);
1506  Rot->rotateY(r_theta * rad);
1507  Rot->rotateZ(rot_z * rad);
1508 
1509  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1510  int copyno = (j_idx << 16) + k_idx;
1511 
1512  new G4PVPlacement(Rot, Crystal_Center,
1513  thirtytwo_logic,
1514  FourByFourName,
1515  ecalenvelope,
1516  0, copyno, OverlapCheck());
1517  }
1518  else if (MappingIndex == 22)
1519  {
1520  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1521 
1522  rot_z = Crystals[j][7] * (-1.0);
1523 
1524  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1525  Rot->rotateX(r_phi * rad);
1526  Rot->rotateY(r_theta * rad);
1527  Rot->rotateZ(rot_z * rad);
1528 
1529  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1530  int copyno = (j_idx << 16) + k_idx;
1531 
1532  new G4PVPlacement(Rot, Crystal_Center,
1533  twentytwo_logic,
1534  FourByFourName,
1535  ecalenvelope,
1536  0, copyno, OverlapCheck());
1537  }
1538  else if (MappingIndex == 12)
1539  {
1540  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1541 
1542  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1543  Rot->rotateX(r_phi * rad);
1544  Rot->rotateY(r_theta * rad);
1545  Rot->rotateZ(rot_z * rad);
1546 
1547  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1548  int copyno = (j_idx << 16) + k_idx;
1549  new G4PVPlacement(Rot, Crystal_Center,
1550  twelve_logic,
1551  FourByFourName,
1552  ecalenvelope,
1553  0, copyno, OverlapCheck());
1554  }
1555 
1556  j++;
1557  }
1558 
1559  //Third Quadrant
1560  j = 0;
1561  while (j_cry > j)
1562  {
1563  MappingIndex = Crystals[j][8];
1564 
1565  j_idx = Crystals[j][0];
1566  k_idx = k_max - Crystals[j][1];
1567  x_cent = Crystals[j][2] - GetParams()->get_double_param("place_x") * cm;
1568  y_cent = -1.0 * (Crystals[j][3] - GetParams()->get_double_param("place_y") * cm);
1569  z_cent = Crystals[j][4] - GetParams()->get_double_param("place_z") * cm;
1570  r_theta = Crystals[j][5];
1571  r_phi = -1.0 * Crystals[j][6];
1572  rot_z = Crystals[j][7] / 4.0;
1573 
1574  if (MappingIndex == 16)
1575  {
1576  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1577 
1578  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1579  Rot->rotateX(r_phi * rad);
1580  Rot->rotateY(r_theta * rad);
1581  Rot->rotateZ(0 * rad);
1582 
1583  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1584  int copyno = (j_idx << 16) + k_idx;
1585 
1586  new G4PVPlacement(Rot, Crystal_Center,
1587  crystal_logic,
1588  FourByFourName,
1589  ecalenvelope,
1590  0, copyno, OverlapCheck());
1591  }
1592  else if (MappingIndex == 32)
1593  {
1594  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1595 
1596  if (j_idx == (k_max - 1) / 2 && k_idx == (k_max - 1) / 2) rot_z = rot_z * 3.0;
1597 
1598  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1599  Rot->rotateX(r_phi * rad);
1600  Rot->rotateY(r_theta * rad);
1601  Rot->rotateZ(rot_z * rad);
1602 
1603  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1604  int copyno = (j_idx << 16) + k_idx;
1605 
1606  new G4PVPlacement(Rot, Crystal_Center,
1607  thirtytwo_logic,
1608  FourByFourName,
1609  ecalenvelope,
1610  0, copyno, OverlapCheck());
1611  }
1612  else if (MappingIndex == 22)
1613  {
1614  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1615 
1616  rot_z = Crystals[j][7] * (1.0);
1617 
1618  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1619  Rot->rotateX(r_phi * rad);
1620  Rot->rotateY(r_theta * rad);
1621  Rot->rotateZ(rot_z * rad);
1622 
1623  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1624  int copyno = (j_idx << 16) + k_idx;
1625  new G4PVPlacement(Rot, Crystal_Center,
1626  twentytwo_logic,
1627  FourByFourName,
1628  ecalenvelope,
1629  0, copyno, OverlapCheck());
1630  }
1631  else if (MappingIndex == 12)
1632  {
1633  G4ThreeVector Crystal_Center = G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1634 
1635  rot_z = rot_z * -4.0000;
1636 
1637  G4RotationMatrix *Rot = new G4RotationMatrix(); //rotation matrix for the placement of each crystal
1638  Rot->rotateX(r_phi * rad);
1639  Rot->rotateY(r_theta * rad);
1640  Rot->rotateZ(rot_z * rad);
1641 
1642  string FourByFourName = "FourByFour_j_" + to_string(j_idx) + "_k_" + to_string(k_idx);
1643  int copyno = (j_idx << 16) + k_idx;
1644 
1645  new G4PVPlacement(Rot, Crystal_Center,
1646  twelve_logic,
1647  FourByFourName,
1648  ecalenvelope,
1649  0, copyno, OverlapCheck());
1650  }
1651 
1652  j++;
1653  }
1654 
1655  return 0;
1656 }