Belle II Software  release-05-01-25
CsiCreator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter, Igal Jaegle, Alexandre Beaulieu *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <beast/csi/geometry/CsiCreator.h>
12 #include <beast/csi/simulation/SensitiveDetector.h>
13 #include <beast/csi/geometry/CsiGeometryPar.h>
14 #include <geometry/Materials.h>
15 #include <geometry/CreatorFactory.h>
16 #include <framework/gearbox/GearDir.h>
17 #include <framework/logging/Logger.h>
18 
19 #include <cmath>
20 #include <boost/format.hpp>
21 #include <boost/foreach.hpp>
22 #include <boost/algorithm/string.hpp>
23 #include <G4AssemblyVolume.hh>
24 #include <G4LogicalVolume.hh>
25 
26 //Shapes
27 #include <G4Trap.hh>
28 #include <G4Box.hh>
29 #include <G4SubtractionSolid.hh>
30 
31 //Visualization Attributes
32 #include <G4VisAttributes.hh>
33 
34 #define PI 3.14159265358979323846
35 
36 using namespace std;
37 using namespace boost;
38 
39 
40 namespace Belle2 {
47  namespace csi {
48 
49  // Register the creator
51  geometry::CreatorFactory<CsiCreator> CsiFactory("CSICreator");
52 
53  // add foil thickness //
54  const double avoidov = 1 + 1E-6;
55 
57 
58  CsiCreator::CsiCreator(): m_sensitive(0)
59  {
60  m_sensitive = new SensitiveDetector();
61  }
62 
63  CsiCreator::~CsiCreator()
64  {
65  if (m_sensitive) delete m_sensitive;
66  }
67 
68  void CsiCreator::create(const GearDir& content, G4LogicalVolume& topVolume, geometry::GeometryTypes type)
69  {
70 
71  //Print the type (not used for now)
72  B2DEBUG(200, "CsI Geometry Type: " << type);
73 
74  // Print list of defined materials
75  /*
76  G4NistManager* nistManager = G4NistManager::Instance();
77  cout << *(G4Material::GetMaterialTable()) << endl;
78  */
79 
80  G4Transform3D BrR = G4RotateZ3D(0.0);
81 
82  int nEnc = content.getNumberNodes("/Enclosures/Enclosure");
83 
84  G4AssemblyVolume* assemblyEnclosures = new G4AssemblyVolume();
85  for (int iEnc = 1; iEnc <= nEnc; iEnc++) {
86  BuildEnclosure(content, assemblyEnclosures, "side", iEnc);
87  //
88  }
89 
90  assemblyEnclosures->MakeImprint(&topVolume, BrR);
91 
92  // Show cell IDs and volume names
93  B2INFO("Positions of the individual CsI crystals");
94  CsiGeometryPar* eclp = CsiGeometryPar::Instance();
95  unsigned int i = 0;
96  for (std::vector<G4VPhysicalVolume*>::iterator it = assemblyEnclosures->GetVolumesIterator();
97  i != assemblyEnclosures->TotalImprintedVolumes();
98  ++it, ++i) {
99 
100  G4VPhysicalVolume* volume = *it;
101  string VolumeName = volume->GetName();
102  if (VolumeName.find("Crystal") < string::npos) {
103  B2INFO("Crystal Number " << eclp->CsiVolNameToCellID(VolumeName) <<
104  " placed at (r[cm],[deg],z[cm]) = (" << setprecision(1) << fixed <<
105  volume->GetTranslation().perp() / CLHEP::cm << "," <<
106  volume->GetTranslation().phi() * 180.0 / PI << "," <<
107  volume->GetTranslation().z() / CLHEP::cm << ")");
108  }
109 
110 
111  }// for all physical volumes in the assembly
112  }// create
113 
114  void CsiCreator::PutCrystal(const GearDir& content,
115  G4AssemblyVolume* assembly,
116  G4Transform3D position,
117  int iEnclosure,
118  int iCry)
119  {
120  if (iCry <= 0) return;
121 
122  GearDir counter(content);
123  double foilthickness = counter.getLength("/Wrapping/Thickness") * CLHEP::cm;
124  G4Material* foilMaterial = geometry::Materials::get(counter.getString("/Wrapping/Material"));
125 
126 
127  int nCry = content.getNumberNodes("/EndCapCrystals/EndCapCrystal");
128  if (iCry > nCry) {
129  B2ERROR("CsiCreator: Crystal index too high");
130  return ;
131  }
132 
133  counter.append((format("/EndCapCrystals/EndCapCrystal[%1%]/") % (iCry)).str());
134  double h1 = counter.getLength("K_h1") * CLHEP::cm;
135  double h2 = counter.getLength("K_h2") * CLHEP::cm;
136  double bl1 = counter.getLength("K_bl1") * CLHEP::cm;
137  double bl2 = counter.getLength("K_bl2") * CLHEP::cm;
138  double tl1 = counter.getLength("K_tl1") * CLHEP::cm;
139  double tl2 = counter.getLength("K_tl2") * CLHEP::cm;
140  double alpha1 = counter.getAngle("K_alpha1");
141  double alpha2 = counter.getAngle("K_alpha2");
142  double halflength = counter.getLength("k_HalfLength") * CLHEP::cm;
143 
144  // Read and create material
145  string strMatCrystal = counter.getString("Material", "Air");
146  G4Material* crystalMaterial = geometry::Materials::get(strMatCrystal);
147 
148 
149  G4VisAttributes* CrystalVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 0.0, 1.0));
150 
151  if (strMatCrystal.compare("CsI") == 0) {
152  CrystalVisAtt->SetColour(18.0 / 256, 230.0 / 256, 3.0 / 256);
153  } else if (strMatCrystal.compare("CsI-Tl") == 0) {
154  CrystalVisAtt->SetColour(0.0, 0.5, 1.0);
155  } else if (strMatCrystal.compare("LYSO") == 0) {
156  CrystalVisAtt->SetColour(0.820, 0.148, 0.1875);
157  }
158 
160  double fwtrapangle1 = atan(2 * h1 / (bl1 - tl1)); // the smaller angle of the trap
161  double fwtrapangle2 = atan(2 * h2 / (bl2 - tl2));
162  double foilh1 = h1 + foilthickness;
163  double foilh2 = h2 + foilthickness;
164  double foiltl1 = tl1 + foilthickness * tan(fwtrapangle1 / 2);
165  double foilbl1 = bl1 + foilthickness / tan(fwtrapangle1 / 2);
166  double foiltl2 = tl2 + foilthickness * tan(fwtrapangle2 / 2);
167  double foilbl2 = foiltl2 + (foilbl1 - foiltl1) * foilh2 / foilh1;
168 
169  double foilhalflength = halflength + foilthickness;
171 
172  string cryLogiVolName = (format("Enclosure_%1%_Crystal_%2%") % iEnclosure % iCry).str();
173  G4Trap* CrystalShape = new G4Trap((format("sCrystal_%1%") % iCry).str().c_str(),
174  halflength , 0 , 0, h1 , bl1, tl1 , alpha1 , h2 , bl2, tl2, alpha2);
175  G4LogicalVolume* Crystal = new G4LogicalVolume(CrystalShape, crystalMaterial,
176  cryLogiVolName.c_str(),
177  0, 0, 0);
178 
179  Crystal->SetVisAttributes(CrystalVisAtt);
180  Crystal->SetSensitiveDetector(m_sensitive);
181 
182  //cout << "CSI volume " << CrystalShape->GetCubicVolume() / CLHEP::cm / CLHEP::cm / CLHEP::cm
183  //<< " density " << crystalMaterial->GetDensity() / CLHEP::g * CLHEP::cm * CLHEP::cm * CLHEP::cm << endl;
184 
186  G4Trap* Foilout = new G4Trap((format("Foilout_%1%") % iCry).str().c_str(),
187  foilhalflength , 0 , 0, foilh1, foilbl1,
188  foiltl1, alpha1 , foilh2, foilbl2,
189  foiltl2, alpha2);
190 
191  G4Trap* Foilin = new G4Trap((format("solidEclCrystal_%1%") % iCry).str().c_str(),
192  halflength * avoidov , 0 , 0, h1 * avoidov ,
193  bl1 * avoidov, tl1 * avoidov , alpha1 , h2 * avoidov,
194  bl2 * avoidov, tl2 * avoidov, alpha2);
195  G4SubtractionSolid* FoilShape = new G4SubtractionSolid((format("sFoil_%1%") % iCry).str().c_str(),
196  Foilout, Foilin);
197 
198  G4LogicalVolume* Foil = new G4LogicalVolume(FoilShape, foilMaterial,
199  (format("Foil_%1%") % iCry).str().c_str(),
200  0, 0, 0);
201 
202  G4VisAttributes* FoilVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1, 0.5));
203  Foil->SetVisAttributes(FoilVisAtt);
204 
205  //Hide the foils for now...
206  Foil->SetVisAttributes(G4VisAttributes::GetInvisible());
207 
208  assembly->AddPlacedVolume(Crystal, position);
209  assembly->AddPlacedVolume(Foil, position);
210 
211  return ;
212 
213  }
214 
215  void CsiCreator::BuildEnclosure(const GearDir& content, G4AssemblyVolume* assembly, string side, int iEnclosure)
216  {
217 
218  string gearPath = "Enclosures/Enclosure";
219  int nEnclosures = content.getNumberNodes(gearPath);
220 
221  if (iEnclosure > nEnclosures) {
222  B2ERROR("Enclosure index too high");
223  return ;
224  }
225 
226  // Build the box (same for all)
227  double width = content.getLength("Enclosures/Width") * CLHEP::cm;
228  double length = content.getLength("Enclosures/Length") * CLHEP::cm;
229  double depth = content.getLength("Enclosures/Depth") * CLHEP::cm;
230  double thk = content.getLength("Enclosures/Thickness") * CLHEP::cm;
231  double fold = content.getLength("Enclosures/Fold") * CLHEP::cm;
232  double lidthk = content.getLength("Enclosures/LidThickness") * CLHEP::cm;
233  double halflength = 15.0 * CLHEP::cm;
234  double zshift = 0.5 * length - thk - halflength; /*< Shift of the box along z-axis to make crystal touch the panel **/
235 
236  string strMatEnclosure = content.getString("Enclosures/Material", "5052-Alloy");
237  G4Material* EnclosureMat = geometry::Materials::get(strMatEnclosure);
238 
239  string strMatEncloLid = content.getString("Enclosures/LidMaterial", "5052-Alloy");
240  G4Material* EncloLidMat = geometry::Materials::get(strMatEncloLid);
241 
242  G4Box* outer = new G4Box("Outer", 0.5 * width, 0.5 * depth, 0.5 * length);
243  G4Box* inner = new G4Box("Inner", 0.5 * width - thk, 0.5 * depth - thk, 0.5 * length - thk);
244  G4Box* opening = new G4Box("Opening", 0.5 * width - fold, 0.5 * depth, 0.5 * length - fold);
245  G4Box* lid = new G4Box("Lid", 0.5 * width, 0.5 * lidthk, 0.5 * length);
246 
247  G4ThreeVector translation(0, thk, 0);
248  G4Translate3D transform(translation);
249  G4SubtractionSolid* enclosureShapeT = new G4SubtractionSolid("EnclosureShapeT", outer, inner);
250  G4SubtractionSolid* enclosureShape = new G4SubtractionSolid("EnclosureShape",
251  enclosureShapeT, opening, transform);
252 
253  //Thread the strings
254  string enclosurePath = (format("/%1%[%2%]") % gearPath % iEnclosure).str();
255  string logiVolName = (format("%1%Enclosure_%2%") % side % iEnclosure).str();
256  string logiLidVolName = (format("%1%EnclosureLid_%2%") % side % iEnclosure).str();
257 
258  // Connect the appropriate Gearbox path
259  GearDir enclosureContent(content);
260  enclosureContent.append(enclosurePath);
261 
262  // Create logical volumes
263  G4LogicalVolume* logiEnclosure = new G4LogicalVolume(enclosureShape, EnclosureMat, logiVolName, 0, 0, 0);
264  G4LogicalVolume* logiEncloLid = new G4LogicalVolume(lid, EncloLidMat, logiLidVolName, 0, 0, 0);
265 
266  // Read position
267  double PosZ = enclosureContent.getLength("PosZ") * CLHEP::cm;
268  double PosR = enclosureContent.getLength("PosR") * CLHEP::cm;
269  double PosT = enclosureContent.getAngle("PosT") ;
270 
271  // Read Orientation
272  double Phi1 = enclosureContent.getAngle("AngPhi1") ;
273  double Theta = enclosureContent.getAngle("AngTheta") ;
274  double Phi2 = enclosureContent.getAngle("AngPhi2") ;
275 
276  //Read position adjustments from nominal
277  double AdjX = enclosureContent.getLength("ShiftX") * CLHEP::cm;
278  double AdjY = enclosureContent.getLength("ShiftY") * CLHEP::cm;
279  double AdjZ = enclosureContent.getLength("ShiftZ") * CLHEP::cm;
280 
281  G4Transform3D zsh = G4Translate3D(0, 0, zshift);
282  G4Transform3D invzsh = G4Translate3D(0, 0, -zshift);
283  G4Transform3D m1 = G4RotateZ3D(Phi1);
284  G4Transform3D m2 = G4RotateY3D(Theta);
285  G4Transform3D m3 = G4RotateZ3D(Phi2);
286  G4Transform3D position = G4Translate3D(PosR * cos(PosT), PosR * sin(PosT), PosZ);
287  G4Transform3D adjust = G4Translate3D(AdjX, AdjY, AdjZ);
288  G4Transform3D lidpos = G4Translate3D(0, 0.5 * (depth + lidthk), 0);
289 
290  G4Transform3D Tr = position * m3 * m2 * m1;
291  G4Transform3D ZshTr = Tr * zsh;
292  G4Transform3D ZshTrAdj = adjust * ZshTr;
293  G4Transform3D LidTr = ZshTr * lidpos;
294  G4Transform3D LidTrAdj = adjust * LidTr;
295 
296  G4VisAttributes* VisAtt = new G4VisAttributes(G4Colour(1.0, 0.5, 0.0, 0.5));
297  logiEnclosure->SetVisAttributes(VisAtt);
298 
299  G4VisAttributes* LidVisAtt = new G4VisAttributes(G4Colour(0.8, 1.0, 0.4, 0.5));
300  logiEncloLid->SetVisAttributes(LidVisAtt);
301  //logiEncloLid->SetVisAttributes(G4VisAttributes::GetInvisible());
302 
303  B2INFO("CsIBox No. " << iEnclosure << " Nominal pos. (mm): " << ZshTr.getTranslation());
304  B2INFO(" Installed pos. (mm): " << ZshTrAdj.getTranslation());
305  B2INFO(" Rotation matrix : " << ZshTrAdj.getRotation());
306  B2INFO(" ");
307 
308  assembly->AddPlacedVolume(logiEnclosure, ZshTrAdj);
309  assembly->AddPlacedVolume(logiEncloLid, LidTrAdj);
310 
311 
312  int nSlots = enclosureContent.getNumberNodes("CrystalInSlot");
313 
314  for (int iSlot = 1; iSlot <= nSlots; iSlot++) {
315  //Thread the strings
316  string slotPath = (format("/Enclosures/Slot[%1%]") % iSlot).str();
317 
318  GearDir slotContent(content);
319  slotContent.append(slotPath);
320 
321  double SlotX = slotContent.getLength("PosX") * CLHEP::cm;
322  double SlotY = slotContent.getLength("PosY") * CLHEP::cm;
323  double SlotZ = slotContent.getLength("PosZ") * CLHEP::cm;
324 
325 
326  G4Transform3D Pos = G4Translate3D(SlotX, SlotY, SlotZ);
327 
328  int CryID = enclosureContent.getInt((format("/CrystalInSlot[%1%]") % iSlot).str());
329 
330  PutCrystal(content, assembly, adjust * Tr * Pos, iEnclosure, CryID);
331  }
332 
333  return;
334 
335  }
336 
337  } // csi namespace
339 } // Belle2 namespace
Belle2::gearbox::Interface::getAngle
double getAngle(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard angle unit.
Definition: Interface.h:301
Belle2::gearbox::Interface::getInt
int getInt(const std::string &path="") const noexcept(false)
Get the parameter path as a int.
Definition: Interface.cc:70
Belle2::csi::SensitiveDetector
Sensitive Detector implementation of the CSI detector.
Definition: SensitiveDetector.h:32
Belle2::csi::CsiGeometryPar
The Class for CSI Geometry Parameters.
Definition: CsiGeometryPar.h:54
Belle2::geometry::Materials::get
static G4Material * get(const std::string &name)
Find given material.
Definition: Materials.h:65
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::GearDir::append
void append(const std::string &path)
Append something to the current path, modifying the GearDir in place.
Definition: GearDir.h:62
Belle2::csi::CsiGeometryPar::CsiVolNameToCellID
int CsiVolNameToCellID(const G4String VolumeName)
Get Cell Id.
Definition: CsiGeometryPar.cc:151
Belle2::gearbox::Interface::getLength
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:261
Belle2::csi::CsiFactory
geometry::CreatorFactory< CsiCreator > CsiFactory("CSICreator")
Creator creates the CSI geometry.
Belle2::csi::avoidov
const double avoidov
foil inside is a little bit lager than crystal to avoid overlap
Definition: CsiCreator.cc:54
Belle2::geometry::GeometryTypes
GeometryTypes
Flag indiciating the type of geometry to be used.
Definition: GeometryManager.h:39
Belle2::GearDir::getNumberNodes
virtual int getNumberNodes(const std::string &path="") const override
Return the number of nodes a given path will expand to.
Definition: GearDir.h:68