6 #include "GeoModelKernel/GeoMaterial.h"
7 #include "GeoModelKernel/GeoBox.h"
8 #include "GeoModelKernel/GeoTube.h"
9 #include "GeoModelKernel/GeoPcon.h"
10 #include "GeoModelKernel/GeoEllipticalTube.h"
11 #include "GeoModelKernel/GeoShapeSubtraction.h"
12 #include "GeoModelKernel/GeoShapeIntersection.h"
13 #include "GeoModelKernel/GeoShapeUnion.h"
14 #include "GeoModelKernel/GeoShapeShift.h"
15 #include "GeoModelKernel/GeoLogVol.h"
16 #include "GeoModelKernel/GeoNameTag.h"
17 #include "GeoModelKernel/GeoPhysVol.h"
18 #include "GeoModelKernel/GeoFullPhysVol.h"
19 #include "GeoModelKernel/GeoDefinitions.h"
20 #include "GeoModelKernel/Units.h"
21 #include "GeoGenericFunctions/AbsFunction.h"
22 #include "GeoGenericFunctions/Variable.h"
23 #include "GeoGenericFunctions/Sin.h"
24 #include "GeoGenericFunctions/Cos.h"
25 #include "CLHEP/Geometry/Point3D.h"
27 #include "GaudiKernel/SystemOfUnits.h"
78 m_properties(
"ForwardRegionProperties")
86 LogStream << MSG::ERROR <<
": Failed to load magnet properties" <<
endmsg;
115 matName =
"std::Vacuum";
116 const GeoMaterial *vacuum = materialManager->
getMaterial(matName);
122 GeoElement *hydrogen =
new GeoElement(
"Hydrogen",
"H",1.0, 1.010);
123 GeoElement *oxygen =
new GeoElement(
"Oxygen",
"O", 8.0, 16.0);
124 water->add(hydrogen,0.11);
125 water->add(oxygen,0.89);
130 const GeoElement*
C = materialManager->
getElement(
"Carbon");
131 const GeoElement*
N = materialManager->
getElement(
"Nitrogen");
132 const GeoElement* Si = materialManager->
getElement(
"Silicon");
133 const GeoElement*
P = materialManager->
getElement(
"Phosphorus");
134 const GeoElement*
S = materialManager->
getElement(
"Sulfur");
135 const GeoElement* Cr = materialManager->
getElement(
"Chromium");
136 const GeoElement* Mn = materialManager->
getElement(
"Manganese");
137 const GeoElement* Fe = materialManager->
getElement(
"Iron");
138 const GeoElement* Ni = materialManager->
getElement(
"Nickel");
139 const GeoElement* Mo = materialManager->
getElement(
"Molybdenum");
140 const GeoElement* Cu = materialManager->
getElement(
"Copper");
143 const GeoElement* Al = materialManager->
getElement(
"Aluminium");
144 const GeoElement*
B = materialManager->
getElement(
"Boron");
145 const GeoElement* O = materialManager->
getElement(
"Oxygen");
150 copper->add(
const_cast<GeoElement*
> (Cu),1.0);
152 m_MapMaterials.insert(std::pair<std::string,GeoMaterial*>(matName,copper));
155 matName =
"Tungsten";
156 const GeoElement*
W = materialManager->
getElement(
"Wolfram");
158 tungsten->add(
const_cast<GeoElement*
> (
W),1.0);
160 m_MapMaterials.insert(std::pair<std::string,GeoMaterial*>(matName,tungsten));
164 matName =
"GlidCopAL15";
167 double aCu, aAl, aO, aB, aTot;
175 glidcop->add(
const_cast<GeoElement*
> (Cu), aCu/aTot);
176 glidcop->add(
const_cast<GeoElement*
> (Al), aAl/aTot);
177 glidcop->add(
const_cast<GeoElement*
> (O), aO/aTot);
178 glidcop->add(
const_cast<GeoElement*
> (
B), aB/aTot);
180 m_MapMaterials.insert(std::pair<std::string,GeoMaterial*>(matName,glidcop));
186 double aC,aN,aSi,aP,aS,aCr,aMn,aFe,aNi,aMo,Atot;
198 Atot=aFe+aC+aMn+aSi+aP+aS+aCr+aMo+aNi+aN;
200 steel->add(
const_cast<GeoElement*
> (Fe),aFe/Atot);
201 steel->add(
const_cast<GeoElement*
> (
C), aC/Atot);
202 steel->add(
const_cast<GeoElement*
> (Mn),aMn/Atot);
203 steel->add(
const_cast<GeoElement*
> (Si),aSi/Atot);
204 steel->add(
const_cast<GeoElement*
> (
P), aP/Atot);
205 steel->add(
const_cast<GeoElement*
> (
S), aS/Atot);
206 steel->add(
const_cast<GeoElement*
> (Cr),aCr/Atot);
207 steel->add(
const_cast<GeoElement*
> (Mo),aMo/Atot);
208 steel->add(
const_cast<GeoElement*
> (Ni),aNi/Atot);
209 steel->add(
const_cast<GeoElement*
> (
N), aN/Atot);
211 m_MapMaterials.insert(std::pair<std::string,GeoMaterial*>(matName,steel));
217 std::ostringstream strs;
232 int name, zStart, zEnd,
type, xAperture, yAperture, xStart, xEnd, yStart, yEnd, tubeThickness;
245 int lDFSize = loadedDataFile.size();
248 for(
int i=0;
i < lDFSize;
i++)
251 HepGeom::Point3D<double> pointMagStart;
252 HepGeom::Point3D<double> pointMagEnd;
256 bool pointsDefined = !(pointMagStart.distance2() == 0 || pointMagEnd.distance2() == 0);
262 double startZ = pointMagStart[2];
263 double endZ = pointMagEnd[2];
264 double startX = pointMagStart[0];
265 double endX = pointMagEnd[0];
266 double rotationAngle = atan2(endX - startX,endZ - startZ);
273 loadedDataFile[
i][xStart] =
num2str(pointMagStart[0]/1000);
274 loadedDataFile[
i][yStart] =
num2str(pointMagStart[1]/1000);
275 loadedDataFile[
i][zStart] =
num2str(pointMagStart[2]/1000);
276 loadedDataFile[
i-1][zEnd] =
num2str(
sgn(pointMagStart[2])*(abs(pointMagStart[2])-dL)/1000);
278 loadedDataFile[
i][xEnd] =
num2str(pointMagEnd[0]/1000);
279 loadedDataFile[
i][yEnd] =
num2str(pointMagEnd[1]/1000);
280 loadedDataFile[
i+1][zStart] =
num2str(
sgn(pointMagEnd[2])*(abs(pointMagEnd[2])+dL)/1000);
281 loadedDataFile[
i][zEnd] =
num2str(pointMagEnd[2]/1000);
286 double x,
y,
z,halfL,rotationAngle,
r,dL,startZ,endZ,startX,endX,startY,endY;
291 for(
int i=0;
i < lDFSize;
i++)
301 x = (startX + endX)/2;
302 y = (startY + endY)/2;
303 z = (startZ + endZ)/2;
307 rotationAngle = atan2(endX - startX,endZ - startZ);
310 halfL = sqrt((endX - startX)*(endX - startX) + (endZ - startZ)*(endZ - startZ))/2;
318 if(loadedDataFile[
i][
name].
find(
"Mag") != std::string::npos)
324 if(
atoi(loadedDataFile[
i][
type].c_str()) == 0){
326 if(loadedDataFile[
i][
name] ==
"VCDBP.7R1.B"){
328 insertCircularElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), trackEnv);
331 insertCircularElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), fwrPhys);
335 if(
atoi(loadedDataFile[
i][
type].c_str()) == 4){
336 if(loadedDataFile[
i][
name] ==
"VCTYF.A4R1.X")
338 else if(loadedDataFile[
i][
name] ==
"TCL.5R1.B1")
347 insertCircularElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), fwrPhys);
353 if(
atoi(loadedDataFile[
i][
type].c_str()) == 1) {
355 insertEllipticalElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), magEnv);
360 if(loadedDataFile[
i][
name].
find(
"Mag") != std::string::npos)
362 if(loadedDataFile[
i][
name] ==
"LQXAA.1R1MagQ1" || loadedDataFile[
i][
name] ==
"LQXAG.3R1MagQ3")
364 if(loadedDataFile[
i][
name] ==
"LQXBA.2R1MagQ2a" || loadedDataFile[
i][
name] ==
"LQXBA.2R1MagQ2b")
370 if(
atoi(loadedDataFile[
i][
type].c_str()) == 2) {
372 insertXRecticircularElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), magEnv);
376 if(
atoi(loadedDataFile[
i][
type].c_str()) == 3) {
378 insertYRecticircularElement(loadedDataFile[
i][
name],
x,
y,
z, rotationAngle,
atof(loadedDataFile[
i][xAperture].c_str()),
atof(loadedDataFile[
i][yAperture].c_str()), halfL, dL,
atof(loadedDataFile[
i][tubeThickness].c_str()), magEnv);
391 LogStream << MSG::INFO <<
"Constructing forward region model" <<
endmsg;
396 std::vector<std::vector<std::string> > loadedDataFileR = this->
loadDataFile((
char*)
"LSS1Rout.csv",11);
397 std::vector<std::vector<std::string> > loadedDataFileL = this->
loadDataFile((
char*)
"LSS1Lout.csv",11);
409 const GeoTube *fwrTubeL =
new GeoTube(0,2*
Gaudi::Units::m,(endZ-startZ)/2);
414 const GeoShapeShift& fwrTube0 = (*fwrTubeL)<<shiftL;
415 const GeoShapeUnion& fwrTube1 = fwrTube0.add((*fwrTubeR)<<shiftR);
423 const GeoShapeSubtraction& fwrTube2 = fwrTube1.subtract((*alfa)<<shiftAlfaL1).subtract((*alfa)<<shiftAlfaL2).subtract((*alfa)<<shiftAlfaR1).subtract((*alfa)<<shiftAlfaR2);
432 const GeoShapeSubtraction& fwrTube3 = fwrTube2.subtract((*afp)<<shiftAfpL1).subtract((*afp)<<shiftAfpR1).subtract((*afp2)<<shiftAfpL2).subtract((*afp2)<<shiftAfpR2);
438 const GeoShapeSubtraction& fwrTube = fwrTube3.subtract((*zdc)<<shiftZdcL1).subtract((*zdc)<<shiftZdcR1);
441 const GeoLogVol *fwrLog =
new GeoLogVol(
"ForwardRegionGeoModel", &fwrTube,
m_MapMaterials[std::string(
"std::Vacuum")]);
442 GeoPhysVol *fwrPhys =
new GeoPhysVol(fwrLog);
443 GeoNameTag *
tag =
new GeoNameTag(
"ForwardRegionGeoModel");
451 LogStream << MSG::INFO <<
"Forward region model succesfully constructed." <<
endmsg;
462 std::vector<std::vector<std::string> > loadedData;
469 LogStream <<
MSG::DEBUG <<
"File " <<
fileName <<
" not found in run directory, trying to load it from DATAPATH" <<
endmsg;
478 std::vector<std::string>
row (
cols);
483 if(
c !=
'@' &&
c !=
'#' &&
c !=
'*' &&
c !=
'$' &&
c !=
'%')
490 loadedData.push_back(
row);
491 file.ignore(1024,
'\n');
494 file.ignore(1024,
'\n');
496 LogStream << MSG::INFO <<
"File " <<
fileName <<
" succesfully loaded." <<
endmsg;