ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4RPDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4RPDetector.cc
1 //____________________________________________________________________________..
2 //
3 // Roman Pot class
4 //
5 // Written by Alexander Bylinkin and Dhevan Gangadharan
6 //
7 // template generated by CreateG4Subsystem.pl
8 //____________________________________________________________________________..
9 
10 #include "EICG4RPDetector.h"
11 
12 #include <phparameter/PHParameters.h>
13 
14 #include <g4main/PHG4Detector.h>
15 
16 #include <Geant4/G4Box.hh>
17 #include <Geant4/G4Color.hh>
18 #include <Geant4/G4LogicalVolume.hh>
19 #include <Geant4/G4Material.hh>
20 #include <Geant4/G4RotationMatrix.hh>
21 #include <Geant4/G4PVPlacement.hh>
22 #include <Geant4/G4SubtractionSolid.hh>
23 #include <Geant4/G4MultiUnion.hh>
24 #include <Geant4/G4UnionSolid.hh>
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4VisAttributes.hh>
27 #include <Geant4/G4TwoVector.hh>
28 #include <Geant4/G4ExtrudedSolid.hh>
29 
30 #include <TSystem.h>
31 
32 #include <cmath>
33 #include <cstdlib>
34 #include <fstream>
35 #include <iostream>
36 #include <sstream>
37 #include <utility>
38 
39 class G4VSolid;
40 class PHCompositeNode;
41 
42 using namespace std;
43 
44 //_______________________________________________________________
46  PHCompositeNode *Node,
48  const std::string &dnam, const int lyr)
49  : PHG4Detector(subsys, Node, dnam)
50  , m_Params(parameters)
51  , m_Layer(lyr)
52 {
53 }
54 
55 //_______________________________________________________________
57 {
58 
59  if( m_ActivePhysicalVolumesMap.find( volume ) != m_ActivePhysicalVolumesMap.end() ) {
60  return 1;
61  }
62 
63  return 0;
64 }
65 
66 //_______________________________________________________________
68 {
69 
70  if( m_VirtualPhysicalVolumesMap.find( volume ) != m_VirtualPhysicalVolumesMap.end() ) {
71  return 1;
72  }
73 
74  return 0;
75 }
76 
77 //_______________________________________________________________
79 {
80  if (IsInDetector(volume))
81  {
82  return 1;
83  }
84 
85  return -1;
86 }
87 
88 //_______________________________________________________________
90 {
91  // Do not forget to multiply the parameters with their respective CLHEP/G4 unit !
92 
93  if (Verbosity() > 1) { std::cout << "Creating Roman Pots" << std::endl; }
94 
95  int layer = m_Params->get_int_param( "layerNumber" );
96  std::string prefix = "Layer" + std::to_string(layer) + "_";
97 
98  double center_X = m_Params->get_double_param( prefix + "pos_x" ) * cm;
99  double center_Y = m_Params->get_double_param( prefix + "pos_y" ) * cm;
100  double center_Z = m_Params->get_double_param( prefix + "pos_z" ) * cm;
101  double overallSize_X = m_Params->get_double_param( prefix + "size_x" ) * cm;
102  double overallSize_Y = m_Params->get_double_param( prefix + "size_y" ) * cm;
103  double enclosureCenter = m_Params->get_double_param( "FFenclosure_center" ) * cm;
104  double rotAngle = m_Params->get_double_param( prefix + "rot_y") * rad;
105 
106  double sensorWidth = m_Params->get_double_param( "Sensor_size" ) * cm;
107  double sensorDepth = m_Params->get_double_param( prefix + "Si_size_z" ) * cm;
108  double CuDepth = m_Params->get_double_param( prefix + "Cu_size_z" ) * cm;
109  double tenSigma_X = m_Params->get_double_param( prefix + "beamHoleHW_x" ) * cm;
110  double tenSigma_Y = m_Params->get_double_param( prefix + "beamHoleHW_y" ) * cm;
111  // Minimum size is set to 100 microns. Zero causes problems with G4ExtrudedSolid below
112  if( tenSigma_X < 0.01 * cm || tenSigma_Y < 0.01 * cm ) {
113  tenSigma_X = 0.01 * cm;
114  tenSigma_Y = 0.01 * cm;
115  }
116 
117  double virtPlaneDepth = 0.001 * cm;
118 
119  // Derived quantities
120  int NsegmentsOffCenter = int( (tenSigma_X - 0.5*sensorWidth) / sensorWidth + 0.5 );
121  int NhorizontalSegments = 4*NsegmentsOffCenter + 2;
122 
123  double x1 = -sensorWidth/2.0;
124  double x2 = sensorWidth/2.0;
125  double y = tenSigma_Y;
126 
127  // Construct polygon xy beam-hole cutout
128  std::vector<G4TwoVector> polygon;
129 
130  if( NhorizontalSegments == 2 ) { // simple case where sensor size is larger than 10sigma
131  // Cutout given by 10sigma
132  polygon.push_back({-tenSigma_X, y});
133  polygon.push_back({tenSigma_X, y});
134  polygon.push_back({tenSigma_X, -y});
135  polygon.push_back({-tenSigma_X, -y});
136  }
137  else { // case where sensors are placed at different y positions to get closer to beam
138 
139  for( int i = 1; i <= NhorizontalSegments; i++ ) {
140 
141  if( i == 1 ) { } // top center
142  else if( i <= 1 + NsegmentsOffCenter ) { // (top) going right
143  // ellipse equation
144  y = tenSigma_Y * sqrt(fabs( 1 - pow( x1/tenSigma_X, 2) ));
145  x2 = x1 + sensorWidth;
146  }
147  else if( i <= 1 + 2*NsegmentsOffCenter ) { // (bottom) going left
148  if( x1 > tenSigma_X ) { y *= -1; }
149  else { y = -tenSigma_Y * sqrt(fabs( 1 - pow( (x1-sensorWidth)/tenSigma_X, 2) ));}
150  x2 = x1 - sensorWidth;
151  }
152  else if( i == 2 + 2*NsegmentsOffCenter ) { // bottom center
153  y = -tenSigma_Y;
154  x2 = x1 - sensorWidth;
155  }
156  else if( i <= 2 + 3*NsegmentsOffCenter ) { // (bottom) going left
157  // ellipse equation
158  y = -tenSigma_Y * sqrt(fabs( 1 - pow( x1/tenSigma_X, 2) ));
159  x2 = x1 - sensorWidth;
160  }
161  else { // (top) going right
162  if( x1 < -tenSigma_X ) { y *= -1; }
163  else { y = tenSigma_Y * sqrt(fabs( 1 - pow( (x1+sensorWidth)/tenSigma_X, 2) )); }
164  x2 = x1 + sensorWidth;
165  }
166 
167  polygon.push_back({x1, y});
168  polygon.push_back({x2, y});
169  x1 = x2;
170  }
171  }
172 
173 
174  // Extrude the xy polygon in z based on sensorDepth (a little more than 1/2 depth)
175  std::vector<G4ExtrudedSolid::ZSection> zsections = {
176  {-0.51*sensorDepth, {0,0}, 1.0}, {+0.51*sensorDepth, {0,0}, 1.0} };
177 
178  G4ExtrudedSolid *polygonCutOut = new G4ExtrudedSolid("Extruded", polygon, zsections);
179 
180  G4Box *FullPlate = new G4Box("FullPlate",overallSize_X/2., overallSize_Y/2., sensorDepth/2.);
181 
182  G4SubtractionSolid *solidRP = new G4SubtractionSolid("EICG4RPSolid", FullPlate, polygonCutOut);
183 
184  G4LogicalVolume *logicalRP = new G4LogicalVolume( solidRP,
185  GetDetectorMaterial( "G4_Si" ), "EICG4RPLogical");
186 
187  G4VisAttributes *vis = new G4VisAttributes( G4Color(1.0, 1.0, 0.0, 1.0) );
188  vis->SetForceSolid(true);
189  logicalRP->SetVisAttributes(vis);
190 
191  G4RotationMatrix *rotm = new G4RotationMatrix();
192  rotm->rotateY( rotAngle );
193 
195  center_X,
196  center_Y,
197  center_Z - enclosureCenter );
198 
199  G4VPhysicalVolume *physicalRP = new G4PVPlacement( rotm, position, logicalRP, "EICG4RP",
200  logicWorld, 0, false, OverlapCheck());
201 
202  // Add it to the list of active volumes so the IsInDetector method picks them up
203  m_ActivePhysicalVolumesMap.insert({physicalRP, layer + 1});
204 
206  // Cu cooling/readout
207 
208  G4Box *CuPlate = (G4Box*)FullPlate->Clone();
209  CuPlate->SetZHalfLength( CuDepth/2.0 );
210  std::vector<G4ExtrudedSolid::ZSection> zsectionsCu = {
211  {-0.51*CuDepth, {0,0}, 1.0}, {+0.51*CuDepth, {0,0}, 1.0} };
212  G4ExtrudedSolid *polygonCutOutCu = new G4ExtrudedSolid("ExtrudedCu", polygon, zsectionsCu);
213  G4SubtractionSolid *solidCu = new G4SubtractionSolid( "EICG4RPCuSolid", CuPlate, polygonCutOutCu);
214  G4LogicalVolume *logicalCu = new G4LogicalVolume( solidCu,
215  GetDetectorMaterial( "G4_Cu" ), "EICG4RPCuLogical");
216 
217  G4VisAttributes *visCu = new G4VisAttributes( G4Color(1.0, 0.0, 1.0, 0.5) );
218  visCu->SetForceSolid(true);
219  logicalCu->SetVisAttributes(visCu);
220  G4ThreeVector positionCu = G4ThreeVector(
221  center_X,
222  center_Y,
223  center_Z + sensorDepth + CuDepth/2.0 - enclosureCenter );
224 
225  G4VPhysicalVolume *physicalCu = new G4PVPlacement( rotm, positionCu, logicalCu, "EICG4RPCu",
226  logicWorld, 0, false, OverlapCheck());
227 
228  // Add it as a passive volume (secondaries created by no hits stored)
229  m_PassivePhysicalVolumesSet.insert(physicalCu);
230 
232  // Virtual plane with no beam hole. In front of active layer
233  G4Box *VirtPlate = new G4Box("VirtPlate",overallSize_X/2., overallSize_Y/2., virtPlaneDepth/2.);
234  G4LogicalVolume *logicalVirt = new G4LogicalVolume( VirtPlate,
235  GetDetectorMaterial( "G4_Galactic" ), "EICG4RPVirtualLogical");
236 
237  G4VisAttributes *visVirt = new G4VisAttributes( G4Color(1.0, 1.0, 1.0, 0.0) );
238  visVirt->SetForceSolid(true);
239  logicalVirt->SetVisAttributes(visVirt);
240 
241  G4ThreeVector positionVirt = G4ThreeVector(
242  center_X,
243  center_Y,
244  center_Z - sensorDepth/2.0 - virtPlaneDepth - enclosureCenter );
245 
246  G4VPhysicalVolume *physicalVirt = new G4PVPlacement( rotm, positionVirt, logicalVirt, "EICG4RPVirt",
247  logicWorld, 0, false, OverlapCheck());
248 
249  // Add it to the virtual volume
250  m_VirtualPhysicalVolumesMap.insert({physicalVirt, layer + 1});
251 
252  return;
253 }
254 
255 //_______________________________________________________________
256 void EICG4RPDetector::Print(const std::string &what) const
257 {
258  std::cout << "EICG4RP Detector:" << std::endl;
259  if (what == "ALL" || what == "VOLUME")
260  {
261  std::cout << "Version 0.1" << std::endl;
262  std::cout << "Parameters:" << std::endl;
263  m_Params->Print();
264  }
265  return;
266 }
267 
269 {
270  return m_Params;
271 }