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>
20 #include <boost/format.hpp>
21 #include <boost/foreach.hpp>
22 #include <boost/algorithm/string.hpp>
23 #include <G4AssemblyVolume.hh>
24 #include <G4LogicalVolume.hh>
29 #include <G4SubtractionSolid.hh>
32 #include <G4VisAttributes.hh>
34 #define PI 3.14159265358979323846
37 using namespace boost;
51 geometry::CreatorFactory<CsiCreator>
CsiFactory(
"CSICreator");
58 CsiCreator::CsiCreator(): m_sensitive(0)
63 CsiCreator::~CsiCreator()
65 if (m_sensitive)
delete m_sensitive;
72 B2DEBUG(200,
"CsI Geometry Type: " << type);
80 G4Transform3D BrR = G4RotateZ3D(0.0);
82 int nEnc = content.getNumberNodes(
"/Enclosures/Enclosure");
84 G4AssemblyVolume* assemblyEnclosures =
new G4AssemblyVolume();
85 for (
int iEnc = 1; iEnc <= nEnc; iEnc++) {
86 BuildEnclosure(content, assemblyEnclosures,
"side", iEnc);
90 assemblyEnclosures->MakeImprint(&topVolume, BrR);
93 B2INFO(
"Positions of the individual CsI crystals");
96 for (std::vector<G4VPhysicalVolume*>::iterator it = assemblyEnclosures->GetVolumesIterator();
97 i != assemblyEnclosures->TotalImprintedVolumes();
100 G4VPhysicalVolume* volume = *it;
101 string VolumeName = volume->GetName();
102 if (VolumeName.find(
"Crystal") < string::npos) {
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 <<
")");
114 void CsiCreator::PutCrystal(
const GearDir& content,
115 G4AssemblyVolume* assembly,
116 G4Transform3D position,
120 if (iCry <= 0)
return;
123 double foilthickness = counter.getLength(
"/Wrapping/Thickness") * CLHEP::cm;
127 int nCry = content.getNumberNodes(
"/EndCapCrystals/EndCapCrystal");
129 B2ERROR(
"CsiCreator: Crystal index too high");
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;
145 string strMatCrystal = counter.getString(
"Material",
"Air");
149 G4VisAttributes* CrystalVisAtt =
new G4VisAttributes(G4Colour(1.0, 1.0, 0.0, 1.0));
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);
160 double fwtrapangle1 = atan(2 * h1 / (bl1 - tl1));
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;
169 double foilhalflength = halflength + foilthickness;
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(),
179 Crystal->SetVisAttributes(CrystalVisAtt);
180 Crystal->SetSensitiveDetector(m_sensitive);
186 G4Trap* Foilout =
new G4Trap((format(
"Foilout_%1%") % iCry).str().c_str(),
187 foilhalflength , 0 , 0, foilh1, foilbl1,
188 foiltl1, alpha1 , foilh2, foilbl2,
191 G4Trap* Foilin =
new G4Trap((format(
"solidEclCrystal_%1%") % iCry).str().c_str(),
195 G4SubtractionSolid* FoilShape =
new G4SubtractionSolid((format(
"sFoil_%1%") % iCry).str().c_str(),
198 G4LogicalVolume* Foil =
new G4LogicalVolume(FoilShape, foilMaterial,
199 (format(
"Foil_%1%") % iCry).str().c_str(),
202 G4VisAttributes* FoilVisAtt =
new G4VisAttributes(G4Colour(0.1, 0.1, 0.1, 0.5));
203 Foil->SetVisAttributes(FoilVisAtt);
206 Foil->SetVisAttributes(G4VisAttributes::GetInvisible());
208 assembly->AddPlacedVolume(Crystal, position);
209 assembly->AddPlacedVolume(Foil, position);
215 void CsiCreator::BuildEnclosure(
const GearDir& content, G4AssemblyVolume* assembly,
string side,
int iEnclosure)
218 string gearPath =
"Enclosures/Enclosure";
219 int nEnclosures = content.getNumberNodes(gearPath);
221 if (iEnclosure > nEnclosures) {
222 B2ERROR(
"Enclosure index too high");
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;
236 string strMatEnclosure = content.getString(
"Enclosures/Material",
"5052-Alloy");
239 string strMatEncloLid = content.getString(
"Enclosures/LidMaterial",
"5052-Alloy");
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);
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);
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();
259 GearDir enclosureContent(content);
260 enclosureContent.
append(enclosurePath);
263 G4LogicalVolume* logiEnclosure =
new G4LogicalVolume(enclosureShape, EnclosureMat, logiVolName, 0, 0, 0);
264 G4LogicalVolume* logiEncloLid =
new G4LogicalVolume(lid, EncloLidMat, logiLidVolName, 0, 0, 0);
267 double PosZ = enclosureContent.
getLength(
"PosZ") * CLHEP::cm;
268 double PosR = enclosureContent.
getLength(
"PosR") * CLHEP::cm;
269 double PosT = enclosureContent.
getAngle(
"PosT") ;
272 double Phi1 = enclosureContent.
getAngle(
"AngPhi1") ;
273 double Theta = enclosureContent.
getAngle(
"AngTheta") ;
274 double Phi2 = enclosureContent.
getAngle(
"AngPhi2") ;
277 double AdjX = enclosureContent.
getLength(
"ShiftX") * CLHEP::cm;
278 double AdjY = enclosureContent.
getLength(
"ShiftY") * CLHEP::cm;
279 double AdjZ = enclosureContent.
getLength(
"ShiftZ") * CLHEP::cm;
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);
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;
296 G4VisAttributes* VisAtt =
new G4VisAttributes(G4Colour(1.0, 0.5, 0.0, 0.5));
297 logiEnclosure->SetVisAttributes(VisAtt);
299 G4VisAttributes* LidVisAtt =
new G4VisAttributes(G4Colour(0.8, 1.0, 0.4, 0.5));
300 logiEncloLid->SetVisAttributes(LidVisAtt);
303 B2INFO(
"CsIBox No. " << iEnclosure <<
" Nominal pos. (mm): " << ZshTr.getTranslation());
304 B2INFO(
" Installed pos. (mm): " << ZshTrAdj.getTranslation());
305 B2INFO(
" Rotation matrix : " << ZshTrAdj.getRotation());
308 assembly->AddPlacedVolume(logiEnclosure, ZshTrAdj);
309 assembly->AddPlacedVolume(logiEncloLid, LidTrAdj);
314 for (
int iSlot = 1; iSlot <= nSlots; iSlot++) {
316 string slotPath = (format(
"/Enclosures/Slot[%1%]") % iSlot).str();
319 slotContent.
append(slotPath);
321 double SlotX = slotContent.
getLength(
"PosX") * CLHEP::cm;
322 double SlotY = slotContent.
getLength(
"PosY") * CLHEP::cm;
323 double SlotZ = slotContent.
getLength(
"PosZ") * CLHEP::cm;
326 G4Transform3D Pos = G4Translate3D(SlotX, SlotY, SlotZ);
328 int CryID = enclosureContent.
getInt((format(
"/CrystalInSlot[%1%]") % iSlot).str());
330 PutCrystal(content, assembly, adjust * Tr * Pos, iEnclosure, CryID);