18 #include "GeoModelKernel/GeoFullPhysVol.h"
20 #include "GeoModelKernel/GeoElement.h"
21 #include "GeoModelKernel/GeoIdentifierTag.h"
22 #include "GeoModelKernel/GeoLogVol.h"
23 #include "GeoModelKernel/GeoMaterial.h"
24 #include "GeoModelKernel/GeoNameTag.h"
25 #include "GeoModelKernel/GeoPhysVol.h"
26 #include "GeoModelKernel/GeoSerialIdentifier.h"
27 #include "GeoModelKernel/GeoSerialTransformer.h"
28 #include "GeoModelKernel/GeoTransform.h"
29 #include "GeoModelKernel/GeoVPhysVol.h"
30 #include "GeoModelKernel/GeoXF.h"
34 #include "GeoModelKernel/GeoPcon.h"
35 #include "GeoModelKernel/GeoTubs.h"
36 #include "GeoModelKernel/GeoCons.h"
37 #include "GeoModelKernel/GeoBox.h"
38 #include "GeoModelKernel/GeoTrap.h"
43 #include "CLHEP/Geometry/Transform3D.h"
44 #include "CLHEP/Vector/Rotation.h"
46 #include "GeoGenericFunctions/Abs.h"
47 #include "GeoGenericFunctions/Sin.h"
48 #include "GeoGenericFunctions/Cos.h"
49 #include "GeoGenericFunctions/Sqrt.h"
50 #include "GeoGenericFunctions/ATan.h"
51 #include "GeoGenericFunctions/Rectangular.h"
52 #include "GeoGenericFunctions/Mod.h"
53 #include "GeoGenericFunctions/Variable.h"
54 #include "GeoGenericFunctions/FixedConstant.h"
56 #include "GaudiKernel/PhysicalConstants.h"
57 #include "GaudiKernel/MsgStream.h"
58 #include "GaudiKernel/Bootstrap.h"
74 #define BUILD_FRONT_ELECTRONICS
75 #define BUILD_HIGHETA_ELECTRONICS
76 #define BUILD_FRONT_G10
77 #define BUILD_BACK_G10
78 #define BUILD_ACCORDION_PLATES
79 #define BUILD_FRONT_STEEL
122 if (!m_ecamPhysicalPos) MakeEnvelope();
123 return m_ecamPhysicalPos;
127 if (!m_ecamPhysicalNeg) MakeEnvelope();
128 return m_ecamPhysicalNeg;
133 ISvcLocator *svcLocator = Gaudi::svcLocator();
135 if (svcLocator->service(
"MessageSvc",
msgSvc,
true )==StatusCode::FAILURE) {
136 throw std::runtime_error(
"Error in BarrelConstruction, cannot access MessageSvc");
138 MsgStream
log(
msgSvc,
"BarrelConstruction");
143 log <<
"++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
144 log <<
"+ +" << std::endl;
145 log <<
"+ Start of Barrel EM GeoModel definition +" << std::endl;
146 log <<
"+ +" << std::endl;
147 log <<
"++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
149 log <<
" " << std::endl;
150 log <<
"Sagging in geometry " << m_A_SAGGING << std::endl;
151 log <<
" " << std::endl;
158 bool doDetailedAbsorberStraight =
false;
159 bool doDetailedAbsorberFold =
false;
163 if(svcLocator->service (
"GeoModelSvc",geoModel) == StatusCode::FAILURE)
164 throw std::runtime_error(
"Error in BarrelConstruction, cannot access GeoModelSvc");
165 if(svcLocator->service (
"RDBAccessSvc",rdbAccess) == StatusCode::FAILURE)
166 throw std::runtime_error(
"Error in BarrelConstruction, cannot access RDBAccessSvc");
169 log << MSG::INFO <<
"Getting primary numbers for " << larVersionKey.
node() <<
", " << larVersionKey.
tag() <<
endmsg;
172 if ((*switchSet).size() !=0) {
175 if (switches->
getInt(
"DETAILED_ABSORBER") !=0) {
177 doDetailedAbsorberStraight =
true;
178 doDetailedAbsorberFold =
true;
186 log << MSG::WARNING <<
" LArSwitches structure not found " <<
endmsg;
189 log << MSG::INFO <<
" Makes detailed absorber sandwich ? " << doDetailedAbsorberStraight <<
" " << doDetailedAbsorberFold <<
endmsg;
190 log << MSG::INFO <<
" Use sagging in geometry ? " << m_A_SAGGING <<
endmsg;
194 GeoGenfun::Sqrt Sqrt;
195 GeoGenfun::ATan ATan;
198 double twopi32 = 2.*twopi64;
202 if (svcLocator->service(
"DetectorStore",
detStore,
false )==StatusCode::FAILURE) {
203 throw std::runtime_error(
"Error in LArDetectorFactory, cannot access DetectorStore");
208 if (StatusCode::SUCCESS !=
detStore->retrieve(materialManager, std::string(
"MATERIALS"))) {
209 throw std::runtime_error(
"Error in BarrelConstruction, stored MaterialManager is not found.");
212 const GeoMaterial *Iron = materialManager->
getMaterial(
"std::Iron");
213 if (!Iron)
throw std::runtime_error(
"Error in BarrelConstruction, std::Iron is not found.");
215 const GeoMaterial *
LAr = materialManager->
getMaterial(
"std::LiquidArgon");
216 if (!
LAr)
throw std::runtime_error(
"Error in BarrelConstruction, std::LiquidArgon is not found.");
218 const GeoMaterial *Lead = materialManager->
getMaterial(
"std::Lead");
219 if (!Lead)
throw std::runtime_error(
"Error in BarrelConstruction, std::Lead is not found.");
221 const GeoMaterial *G10_bar = materialManager->
getMaterial(
"LAr::G10_bar");
222 if (!G10_bar)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::G10_bar is not found.");
224 const GeoMaterial *Moth_elect = materialManager->
getMaterial(
"LAr::MBoards");
225 if (!Moth_elect)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::MBoards is not found.");
227 const GeoMaterial *Cable_elect = materialManager->
getMaterial(
"LAr::Cables");
228 if (!Cable_elect)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::Cables is not found.");
230 const GeoMaterial *Thin_abs = materialManager->
getMaterial(
"LAr::Thinabs");
231 if (!Thin_abs)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::Thinabs is not found.");
233 const GeoMaterial *Thick_abs = materialManager->
getMaterial(
"LAr::Thickabs");
234 if (!Thick_abs)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::Thickabs is not found.");
236 const GeoMaterial *Kapton_Cu = materialManager->
getMaterial(
"LAr::KaptonC");
237 if (!Kapton_Cu)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::KaptonC is not found.");
239 const GeoMaterial *Sumb = materialManager->
getMaterial(
"LAr::SBoard");
240 if (!Sumb)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::SBoard is not found.");
242 const GeoMaterial *Glue = materialManager->
getMaterial(
"LAr::Glue");
243 if (!Glue)
throw std::runtime_error(
"Error in BarrelConstruction, LAr::Glue is not found.");
245 const GeoElement *Pb = materialManager->
getElement(
"Lead");
246 if (!Pb)
throw std::runtime_error(
"Error in BarrelConstruction, element lead not found ");
248 const GeoElement *Fe = materialManager->
getElement(
"Iron");
249 if (!Fe)
throw std::runtime_error(
"Error in BarrelConstruction, element Fe not found ");
252 double contract = m_parameters->GetValue(
"LArEMBAbsorberContraction");
254 double density_lead = Lead->getDensity()/(contract*contract*contract);
259 double density_iron = Iron->getDensity()/(contract*contract*contract);
267 std::cout <<
"*** BarrelConstruction: List of Materials, density,X0,Lambda " << std::endl;
268 std::cout <<
" Iron " << Iron->getDensity() <<
" "
269 << Iron->getRadLength() <<
" "
270 << Iron->getIntLength() << std::endl;
271 std::cout <<
" LAR " <<
LAr->getDensity() <<
" "
272 <<
LAr->getRadLength() <<
" "
273 <<
LAr->getIntLength() << std::endl;
274 std::cout <<
" G10_bar " << G10_bar->getDensity() <<
" "
275 << G10_bar->getRadLength() <<
" "
276 << G10_bar->getIntLength() << std::endl,
277 std::cout <<
" Moth_elect " << Moth_elect->getDensity() <<
" "
278 << Moth_elect->getRadLength() <<
" "
279 << Moth_elect->getIntLength() << std::endl;
280 std::cout <<
" Cable " << Cable_elect->getDensity() <<
" "
281 << Cable_elect->getRadLength() <<
" "
282 << Cable_elect->getIntLength() << std::endl;
283 std::cout <<
" Sumb " << Sumb->getDensity() <<
" "
284 << Sumb->getRadLength() <<
" "
285 << Sumb->getIntLength() << std::endl;
286 std::cout <<
" Thin_abs " << Thin_abs->getDensity() <<
" "
287 << Thin_abs->getRadLength() <<
" "
288 << Thin_abs->getIntLength() << std::endl;
289 std::cout <<
" Thick_abs " << Thick_abs->getDensity() <<
" "
290 << Thick_abs->getRadLength() <<
" "
291 << Thick_abs->getIntLength() << std::endl;
292 std::cout <<
" Kapton_Cu " << Kapton_Cu->getDensity() <<
" "
293 << Kapton_Cu->getRadLength() <<
" "
294 << Kapton_Cu->getIntLength() << std::endl;
296 std::string baseName =
"LAr::EMB::";
303 double Moth_Z_min = m_parameters->GetValue(
"LArEMBMotherZmin");
304 double Moth_Z_max = m_parameters->GetValue(
"LArEMBMotherZmax");
308 double Moth_inner_radius = m_parameters->GetValue(
"LArEMBMotherRmin");
309 double Moth_outer_radius = m_parameters->GetValue(
"LArEMBMotherRmax");
311 double Moth_Phi_Min = 0.;
312 double Moth_Phi_Max = m_parameters->GetValue(
"LArEMBphiMaxBarrel")*
Gaudi::Units::deg;
315 std::cout <<
" *** Mother volume (Ecam) parameters " << std::endl;
316 std::cout <<
" Rmin/Rmax " << Moth_inner_radius <<
" " << Moth_outer_radius << std::endl;
317 std::cout <<
" Zmin/Zmax " << Moth_Z_min <<
" " << Moth_Z_max << std::endl;
322 int Nbrt = (
int) ( m_parameters->GetValue(
"LArEMBnoOFAccZigs") );
325 double Bar_Z_min = Moth_Z_min;
326 double Bar_Z_max = Moth_Z_max;
330 double Rhocen1 = m_parameters->GetValue(
"LArEMBRadiusAtCurvature",0);
333 double Bar_Eta_cut = (
double) (m_parameters->GetValue(
"LArEMBMaxEtaAcceptance"));
335 std::cout <<
" Bar_Eta_cut " << Bar_Eta_cut << std::endl;
338 double Bar_Z_cut, Bar_Rcmx ;
340 Bar_Z_cut = (
double) (m_parameters->GetValue(
"LArEMBMotherZmin"))
341 + (
double) (m_parameters->GetValue(
"LArEMBG10FrontDeltaZ"));
343 Bar_Rcmx = (
double) (m_parameters->GetValue(
"LArEMBRminHighZ"));
348 double Zhalf = (Bar_Z_max-Bar_Z_min)/2.;
350 double Zhalfc = (Bar_Z_cut-Bar_Z_min)/2.;
358 double Rhocen[15], Phicen[15], Delta[15], deltay[15], deltax[15], TetaTrans[15];
359 for (
int idat=0; idat<15; idat++)
362 (
double) (m_parameters->GetValue(
"LArEMBRadiusAtCurvature",idat));
364 (
double) (m_parameters->GetValue(
"LArEMBPhiAtCurvature",idat));
366 (
double) (m_parameters->GetValue(
"LArEMBDeltaZigAngle",idat));
374 (
double) (m_parameters->GetValue(
"LArEMBSaggingAmplitude",idat));
376 (
double) (m_parameters->GetValue(
"LArEMBSaggingAmplitude2",idat));
384 double etaTrans = (
double) (m_parameters->GetValue(
"LArEMBEtaTrans",idat));
385 TetaTrans[idat] = 2.*
atan(
exp(-etaTrans));
388 log <<
MSG::DEBUG <<
"idat " << idat <<
" Rhocen/Phice/Delta/deltay/deltax/etatrans "
390 << Delta[idat]*(1./
Gaudi::Units::deg) <<
" " << deltay[idat] <<
" " << deltax[idat]
391 <<
" " << etaTrans <<
endmsg;
399 if (Phicen[0]<0.) checkParity=1;
402 int Ncell = (
int) (m_parameters->GetValue(
"LArEMBnoOFPhysPhiCell"));
403 int Nabsorber = (Ncell==64) ? Ncell + 1 : Ncell;
404 int Nelectrode = Ncell;
408 if (m_NVISLIM > 0) Nabsorber = m_NVISLIM;
409 if (m_NVISLIM > 0) Nelectrode = m_NVISLIM;
412 log <<
"================LAr BarrelConstruction===========================" << std::endl;
413 log <<
"===YOU ARE RUNNING WITH A LIMITED SLICE OF BARREL ACCORDEON ===" << std::endl;
414 log <<
"===THIS OPTION IS FOR VISUALIZATION OR OTHER TESTING, ONLY!! ===" << std::endl;
415 log <<
"===AND IF THIS IS NOT YOUR INTENTION PLEASE BE SURE TO SET ===" << std::endl;
416 log <<
"===THE FLAG GeoModelSvc.LArDetectorTool.BarrelCellVisLimit=-1 ===" << std::endl;
417 log <<
"===Remember there are no defaults in Athena. Thank you. ===" << std::endl;
418 log <<
"=================================================================" << std::endl;
424 double z1 = m_parameters->GetValue(
"LArEMBG10FrontDeltaZ")+ Moth_Z_min;
426 const double r1 = Moth_inner_radius
427 + m_parameters->GetValue(
"LArEMBInnerElectronics")
428 + m_parameters->GetValue(
"LArEMBG10SupportBarsIn");
430 const double r2 = m_parameters->GetValue(
"LArEMBRminHighZ");
431 const double inv_r2_r1 = 1. / (r2-r1);
434 double zmax1_Stac = z1 + (Rhocen[0]-r1)*(Moth_Z_max-z1)*inv_r2_r1;
452 double Zplan[] = {Bar_Z_min-Zp0,Bar_Z_cut-Zp0,Bar_Z_max-Zp0};
453 double Riacc[] = {Moth_inner_radius,Moth_inner_radius, Rhocen1-safety_rhocen1};
454 double Roacc[] = {Moth_outer_radius,Moth_outer_radius,Moth_outer_radius};
457 std::cout <<
" *** Parameters for ECAM Pcon volume " << std::endl;
458 std::cout <<
" Zplan " << Zplan[0] <<
" " << Zplan[1] <<
" " << Zplan[2] << std::endl;
459 std::cout <<
" Rin " << Riacc[0] <<
" " << Riacc[1] <<
" " << Riacc[2] << std::endl;
460 std::cout <<
" Rout " << Roacc[0] <<
" " << Roacc[1] <<
" " << Roacc[2] << std::endl;
463 int ecamArraySize =
sizeof(Zplan) /
sizeof(
double);
464 std::string
name = baseName +
"ECAM";
465 GeoPcon* pCon =
new GeoPcon(Moth_Phi_Min2,
468 for (
int i=0;
i< ecamArraySize;
i++)
470 pCon->addPlane(Zplan[
i],Riacc[
i],Roacc[
i]);
472 const GeoLogVol* logVol =
new GeoLogVol(
name,pCon,
LAr);
473 m_ecamPhysicalPos =
new GeoFullPhysVol(logVol);
474 m_ecamPhysicalNeg =
new GeoFullPhysVol(logVol);
482 if(!
status.isSuccess())
throw std::runtime_error (
"Cannot store EMB_POS");
487 if(!
status.isSuccess())
throw std::runtime_error (
"Cannot store EMB_NEG");
492 if(!m_fullGeo)
return;
497 #ifdef BUILD_FRONT_ELECTRONICS
503 GeoIntrusivePtr<GeoPhysVol>Elnicsf_phys=
nullptr;
508 Xel1f = m_parameters->GetValue(
"LArEMBInnerElectronics");
509 double DeltaZ = Zhalfc;
510 double Zpos = Zhalfc+Bar_Z_min;
513 std::string
name = baseName +
"TELF";
515 std::cout <<
" *** parameters for TELF tubs " << std::endl;
516 std::cout <<
" DeltaZ " << DeltaZ << std::endl;
517 std::cout <<
" Rmin/Rmax " << Rmini <<
" " << Rmaxi << std::endl,
519 std::cout <<
" Zpos in ECAM " << Zpos << std::endl;
521 GeoTubs* tubs =
new GeoTubs(Rmini,
526 const GeoLogVol* logVol =
new GeoLogVol(
name,tubs,
LAr);
527 Elnicsf_phys =
new GeoPhysVol(logVol);
528 m_ecamPhysicalPos->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
529 m_ecamPhysicalPos->add(Elnicsf_phys);
530 m_ecamPhysicalNeg->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
531 m_ecamPhysicalNeg->add(Elnicsf_phys);
539 GeoIntrusivePtr<GeoPhysVol>Sumb_phys=
nullptr;
542 double Rmini = Moth_inner_radius+Xel1f-ThickSum;
544 double DeltaZ = Zhalfc;
546 std::string
name = baseName +
"SUMB";
548 std::cout <<
" *** parameters for SUMB tubs " << std::endl;
549 std::cout <<
" DeltaZ " << DeltaZ << std::endl;
550 std::cout <<
" Rmin/Rmax " << Rmini <<
" " << Rmaxi << std::endl,
552 std::cout <<
" Zpos in TELF " << Zpos << std::endl;
555 GeoTubs * tubs =
new GeoTubs(Rmini,Rmaxi,DeltaZ,Moth_Phi_Min,Moth_Phi_Max);
556 const GeoLogVol* logVol =
new GeoLogVol(
name,tubs,Sumb);
557 Sumb_phys =
new GeoPhysVol(logVol);
558 Elnicsf_phys->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
559 Elnicsf_phys->add(Sumb_phys);
565 double ClearancePS = m_parameters->GetValue(
"LArEMBMoBoclearfrPS");
566 double RhoPosB = Moth_inner_radius + ClearancePS;
567 double bdx = .5*(m_parameters->GetValue(
"LArEMBMoBoTchickness"));
568 double bdy = .5*(m_parameters->GetValue(
"LArEMBMoBoHeight"));
575 int NoOFboard = (
int) m_parameters->GetValue(
"LArEMBnoOFmothboard");
576 double PhiPos0 = twopi64 + PhiOrig;
578 std::string
name = baseName +
"MOAC";
580 std::cout <<
" *** parameters for MotherBoard (box)" << std::endl;
581 std::cout <<
" dx,dy,dz " << bdx <<
" " << bdy <<
" " << bdz << std::endl;
582 std::cout <<
" Radius pos " << RhoPosB << std::endl;
585 GeoBox * box =
new GeoBox(bdx,bdy,bdz);
586 const GeoLogVol * logVol =
new GeoLogVol(
name,box,Moth_elect);
587 GeoPhysVol * pV =
new GeoPhysVol(logVol);
588 GeoSerialIdentifier * iD =
new GeoSerialIdentifier(0);
589 GeoGenfun::Variable
c;
590 GeoGenfun::GENFUNCTION PhiPos = PhiPos0 + twopi32*
c;
591 GeoXF::TRANSFUNCTION TX =
592 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),RhoPosB*Cos(PhiPos))*
593 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),RhoPosB*Sin(PhiPos))*
594 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),PhiPos);
595 GeoSerialTransformer *st =
new GeoSerialTransformer(pV, &TX, NoOFboard);
596 Elnicsf_phys->add(iD);
597 Elnicsf_phys->add(st);
613 double Dx1 = .5*(m_parameters->GetValue(
"LArEMBCablethickat0"));
614 double Dx2 = .5*Bar_Eta_cut*(m_parameters->GetValue(
"LArEMBthickincrfac"));
616 double Dy1 = .5*(m_parameters->GetValue(
"LArEMBCableEtaheight"));
619 int NoOFcable = (
int) m_parameters->GetValue(
"LArEMBnoOFcableBundle");
620 double RhoPosC = Moth_inner_radius + ClearancePS;
622 std::string
name = baseName +
"CAAC";
624 std::cout <<
" *** Parameters for Cable (Trap) " << std::endl;
625 std::cout <<
" Dzc " << Dzc << std::endl;
626 std::cout <<
" Dx1,Dx2 (half thickness at eta=0 and etamax " << Dx1 <<
" " << Dx2 << std::endl;
627 std::cout <<
" Dy1,Dy2 " << Dy1 <<
" " << Dy2 << std::endl;
628 std::cout <<
" Radial Position " << RhoPosC << std::endl;
631 GeoTrap* trap =
new GeoTrap(Dzc,0.,0.,Dy1,Dx1,Dx1,0.,
633 const GeoLogVol* logVol =
new GeoLogVol(
name,trap,Cable_elect);
635 GeoPhysVol * pV =
new GeoPhysVol(logVol);
636 GeoSerialIdentifier * iD =
new GeoSerialIdentifier(0);
637 double PhiPos0 = (twopi64 - asin(bdy/RhoPosB))/2.;
639 std::cout <<
" PhiPos0 " << PhiPos0 << std::endl;
641 GeoGenfun::Variable
I;
644 GeoGenfun::Rectangular KDelta;
645 KDelta.baseline().setValue(0.0);
646 KDelta.height().setValue(1.0);
647 KDelta.x0().setValue(-0.5);
648 KDelta.x1().setValue(0.5);
651 GeoGenfun::Mod Mod1(1.0),Mod2(2.0),Mod4(4.0);
652 GeoGenfun::GENFUNCTION Int =
I - Mod1;
653 GeoGenfun::GENFUNCTION Ccopy = Int(
I + 0.5);
655 GeoGenfun::GENFUNCTION PhiPos1 = PhiPos0 + PhiOrig;
656 GeoGenfun::GENFUNCTION PhiPos2 = twopi32 - PhiPos0 + PhiOrig;
657 GeoGenfun::GENFUNCTION PhiPos00 =
658 (KDelta(Mod4(Ccopy)-2) + KDelta(Mod4(Ccopy)-3))*PhiPos2 +
659 (1.0-KDelta(Mod4(Ccopy)-2)-KDelta(Mod4(Ccopy)-3))*PhiPos1;
660 GeoGenfun::GENFUNCTION PhiPos = PhiPos00 + Mod2(Ccopy)*twopi32;
661 GeoXF::TRANSFUNCTION TX =
662 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),RhoPosC*Cos(PhiPos))*
663 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),RhoPosC*Sin(PhiPos))*
664 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),PhiPos);
665 GeoSerialTransformer *st =
new GeoSerialTransformer(pV, &TX, NoOFcable);
667 for (
int ii=0;ii<NoOFcable;ii++) {
668 std::cout <<
"copy, phi " << ii <<
" " << PhiPos(ii)/
Gaudi::Units::deg << std::endl;
671 Elnicsf_phys->add(iD);
672 Elnicsf_phys->add(st);
677 #endif // BUILD_FRONT_ELECTRONICS
684 #ifdef BUILD_HIGHETA_ELECTRONICS
696 double ze1= zmax1_Stac+clearance;
697 double ze2 = Bar_Z_max;
699 double DeltaZ = 0.5*(ze2-ze1)-safety;
700 double Zpos = ze1+DeltaZ+0.5*safety;
705 std::string
name = baseName +
"TELB";
707 std::cout <<
" *** Parameters for high eta electronics (Cons) " <<std::endl;
708 std::cout <<
" Rmini1,Rmini2,Rmaxi1,Rmaxi2 " << Rmini1 <<
" " << Rmini2 <<
" "
709 << Rmaxi1 <<
" " << Rmaxi2 << std::endl,
710 std::cout <<
" DeltaZ " << DeltaZ << std::endl;
712 std::cout <<
" Zpos " << Zpos << std::endl;
714 GeoCons* cons =
new GeoCons(Rmini1,Rmini2,Rmaxi1,Rmaxi2,
715 DeltaZ,Moth_Phi_Min,Moth_Phi_Max);
716 const GeoLogVol* logVol =
new GeoLogVol(
name,cons,Cable_elect);
717 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
718 m_ecamPhysicalPos->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
719 m_ecamPhysicalPos->add(physVol);
720 m_ecamPhysicalNeg->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
721 m_ecamPhysicalNeg->add(physVol);
726 #endif // BUILD_HIGHETA_ELECTRONICS
733 #ifdef BUILD_FRONT_G10
738 double Xel1f = m_parameters->GetValue(
"LArEMBInnerElectronics");
739 double Xg10f = m_parameters->GetValue(
"LArEMBG10SupportBarsIn");
740 double DeltaZ = 0.5* m_parameters->GetValue(
"LArEMBG10FrontDeltaZ");
741 double Zpos = DeltaZ+Bar_Z_min;
742 double Rmini = Moth_inner_radius + Xel1f;
743 double Rmaxi = Rmini+Xg10f;
744 std::string
name = baseName +
"GTENF";
746 std::cout <<
" *** parameters for front G10 ring (tubs) " << std::endl;
747 std::cout <<
" Rmini,Rmaxi " << Rmini <<
" " << Rmaxi << std::endl;
748 std::cout <<
" DeltaZ " << DeltaZ << std::endl;
750 std::cout <<
" Zpos " << Zpos << std::endl;
752 GeoTubs* tubs =
new GeoTubs(Rmini,Rmaxi,DeltaZ,Moth_Phi_Min,Moth_Phi_Max);
753 const GeoLogVol* logVol =
new GeoLogVol(
name,tubs,G10_bar);
754 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
755 m_ecamPhysicalPos->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
756 m_ecamPhysicalPos->add(physVol);
757 m_ecamPhysicalNeg->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
758 m_ecamPhysicalNeg->add(physVol);
762 larVersionKey.
node());
763 if (extraCones->size() > 0 ) {
765 const std::string& conName = cone->getString(
"CONE");
766 if (conName==
"ExtraInBar") {
767 double extra_dz = 0.5*( cone->getDouble(
"DZ") );
768 double extra_rmin1 = cone->getDouble(
"RMIN1");
769 double extra_rmin2 = cone->getDouble(
"RMIN2");
770 double extra_rmax1 = cone->getDouble(
"RMAX1");
771 double extra_rmax2 = cone->getDouble(
"RMAX2");
772 double extra_phi0 = cone->getDouble(
"PHI0");
773 double extra_dphi = cone->getDouble(
"DPHI");
774 double extra_zpos = cone->getDouble(
"ZPOS");
776 std::string
name = baseName +
"ExtraMat1";
777 GeoCons* cons =
new GeoCons(extra_rmin1,extra_rmin2,extra_rmax1,extra_rmax2,
778 extra_dz,extra_phi0,extra_dphi);
779 const GeoLogVol* logVol =
new GeoLogVol(
name,cons,Lead);
780 GeoIntrusivePtr<GeoPhysVol> physVol2 =
new GeoPhysVol(logVol);
781 physVol->add(
new GeoTransform(GeoTrf::TranslateZ3D(extra_zpos)));
782 physVol->add(physVol2);
789 #endif // BUILD_FRONT_G10
791 #ifdef BUILD_BACK_G10
795 double DeltaZ = Zhalf;
796 double Zpos = Zhalf+Bar_Z_min;
797 double Xtal = m_parameters->GetValue(
"LArEMBLArGapTail")+ 0.1*
Gaudi::Units::mm;
798 double Rmini = Rhocen[Nbrt]+Xtal;
803 std::string
name = baseName +
"GTENB";
805 std::cout <<
" *** parameters for back G10 ring (tubs) " << std::endl;
806 std::cout <<
" Rmini,Rmaxi " << Rmini <<
" " << Rmaxi << std::endl;
807 std::cout <<
" DeltaZ " << DeltaZ << std::endl;
809 std::cout <<
" Zpos " << Zpos << std::endl;
812 GeoTubs* tubs =
new GeoTubs(Rmini,Rmaxi,DeltaZ,Moth_Phi_Min,Moth_Phi_Max);
813 const GeoLogVol* logVol =
new GeoLogVol(
name,tubs,G10_bar);
814 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
815 m_ecamPhysicalPos->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
816 m_ecamPhysicalPos->add(physVol);
817 m_ecamPhysicalNeg->add(
new GeoTransform(GeoTrf::TranslateZ3D(Zpos)));
818 m_ecamPhysicalNeg->add(physVol);
821 #endif // BUILD_BACK_G10
824 GeoIntrusivePtr<GeoPhysVol>stacPhysical=
nullptr;
835 double Zplan1[] = {Bar_Z_min,zmax1_Stac+clearance,Bar_Z_max};
836 double Riacc1[] = {Rhocen[0],Rhocen[0], Bar_Rcmx-clearance};
837 double Roacc1[] = {Rhocen[Nbrt],Rhocen[Nbrt],Rhocen[Nbrt]};
839 std::string
name = baseName +
"STAC";
841 std::cout <<
" *** Parameters for STAC Pcon volume " << std::endl;
842 std::cout <<
" Zplan " << Zplan1[0] <<
" " << Zplan1[1] <<
" " << Zplan1[2] << std::endl;
843 std::cout <<
" Rin " << Riacc1[0] <<
" " << Riacc1[1] <<
" " << Riacc1[2] << std::endl;
844 std::cout <<
" Rout " << Roacc1[0] <<
" " << Roacc1[1] <<
" " << Roacc1[2] << std::endl;
846 std::cout <<
" Zpos " << -Zp0 << std::endl;
848 GeoPcon* pCone =
new GeoPcon(Moth_Phi_Min2,Moth_Phi_Max2);
849 for (
int i=0;
i<3;
i++) pCone->addPlane(Zplan1[
i],Riacc1[
i],Roacc1[
i]);
850 const GeoLogVol* logVol =
new GeoLogVol(
name,pCone,
LAr);
851 stacPhysical =
new GeoPhysVol(logVol);
852 m_ecamPhysicalPos->add(
new GeoTransform(GeoTrf::TranslateZ3D(-Zp0)));
853 m_ecamPhysicalPos->add(stacPhysical);
854 m_ecamPhysicalNeg->add(
new GeoTransform(GeoTrf::TranslateZ3D(-Zp0)));
855 m_ecamPhysicalNeg->add(stacPhysical);
859 #ifdef BUILD_ACCORDION_PLATES
863 double Xtip_pb = m_parameters->GetValue(
"LArEMBLeadTipThickFront");
864 double Xtip_pc = m_parameters->GetValue(
"LArEMBLeadTipThickEnd");
865 double Xtip_gt = m_parameters->GetValue(
"LArEMBG10TipThickFront");
866 double Xtip_gs = m_parameters->GetValue(
"LArEMBG10TipThickEnd");
867 double Thgl = m_parameters->GetValue(
"LArEMBThickAbsGlue");
868 double Thfe = m_parameters->GetValue(
"LArEMBThickAbsIron");
869 double Thpb = m_parameters->GetValue(
"LArEMBThickAbsLead");
870 double Thpb_thin = m_parameters->GetValue(
"LArEMBThinAbsLead");
871 double Thcu = m_parameters->GetValue(
"LArEMBThickElecCopper");
872 double Thfg = m_parameters->GetValue(
"LArEMBThickElecKapton");
873 double Psi = m_parameters->GetValue(
"LArEMBPhiGapAperture");
874 double Contract = m_parameters->GetValue(
"LArEMBAbsorberContraction");
876 double Thce = (Thpb+Thgl+Thfe)*Contract;
877 double Thel = Thcu+Thfg;
881 double Zmin = Bar_Z_min;
882 double Zmax = Bar_Z_max;
896 double Zcp1l[14]={},Zcp1h[14]={},Zcp2l[14]={},Zcp2h[14]={};
897 double Rhol[14]={},Rhoh[14]={};
904 double Cenx[15]={0}, Ceny[15]={0} ;
905 for (
int jf=0; jf<(Nbrt+1); jf++)
907 Cenx[jf] = Rhocen[jf]*
cos(Phicen[jf]);
908 Ceny[jf] = Rhocen[jf]*
sin(Phicen[jf]);
909 Zcp1[jf] = Rhocen[jf]/
tan(TetaTrans[jf]);
911 Zcp2[jf] = z1 + (Rhocen[jf]-r1)*inv_r2_r1*(Moth_Z_max-z1);
916 std::cout <<
" jf " << jf
917 <<
"Cenx/Ceny " << Cenx[jf] <<
" " << Ceny[jf] <<
" "
918 <<
"Zcp1/Zcp2 " << Zcp1[jf] <<
" " << Zcp2[jf] << std::endl;
924 double Rint = m_parameters->GetValue(
"LArEMBNeutFiberRadius");
925 for (
int jl=0; jl<Nbrt; jl++)
927 double Adisc2 = (Cenx[jl+1]-Cenx[jl])*(Cenx[jl+1]-Cenx[jl])+
928 (Ceny[jl+1]-Ceny[jl])*(Ceny[jl+1]-Ceny[jl]);
929 double Along2 = Adisc2 - 4.*Rint*Rint;
930 Along[jl] = sqrt(Along2);
932 std::cout <<
"jl " << jl
933 <<
"straight length " << Along[jl] << std::endl;
935 double ddelta =
M_PI/2.-Delta[jl];
937 if (checkParity==0) {
938 if (jl%2==1) ddelta=-1.*ddelta;
940 if (jl%2==0) ddelta=-1.*ddelta;
942 double xlong=Along[jl]-2.*safety_along;
943 double x2=0.5*(Cenx[jl+1]+Cenx[jl])+0.5*xlong*
cos(ddelta);
944 double y2=0.5*(Ceny[jl+1]+Ceny[jl])+0.5*xlong*
sin(ddelta);
945 double x1=0.5*(Cenx[jl+1]+Cenx[jl])-0.5*xlong*
cos(ddelta);
946 double y1=0.5*(Ceny[jl+1]+Ceny[jl])-0.5*xlong*
sin(ddelta);
949 Zcp1l[jl] = Rhol[jl]/
tan(TetaTrans[jl]);
950 Zcp1h[jl] = Rhoh[jl]/
tan(TetaTrans[jl+1]);
952 Zcp2l[jl] = z1 + (Rhol[jl]-r1)*inv_r2_r1*(Moth_Z_max-z1);
954 Zcp2l[jl]=Moth_Z_max;
956 Zcp2h[jl] = z1 + (Rhoh[jl]-r1)*inv_r2_r1*(Moth_Z_max-z1);
958 Zcp2h[jl]=Moth_Z_max;
960 std::cout <<
"x1,y1,x2,y2 " <<
x1 <<
" " <<
y1 <<
" " <<
x2 <<
" "
962 std::cout <<
" lower/upper radius " << Rhol[jl] <<
" " << Rhoh[jl]
963 <<
" Z for eta0.8 " << Zcp1l[jl] <<
" " << Zcp1h[jl]
964 <<
" Z for etamax " << Zcp2l[jl] <<
" " << Zcp2h[jl] << std::endl;
970 double Gama0 = m_parameters->GetValue(
"LArEMBAbsPhiFirst");
972 GeoGenfun::Variable icopy;
973 GeoGenfun::GENFUNCTION Game = Gama0 + Psi/2 + Alfa*icopy;
974 GeoGenfun::GENFUNCTION Gama = Gama0 + Alfa*icopy;
979 enum FB {FRONT, BACK, TERM};
981 for (
int fb=FRONT; fb!=TERM; fb++)
984 int irl = fb==FRONT ? 0 : Nbrt;
986 double Xtips = Xtip_pb + Xtip_gt;
987 double Xtipt = Xtip_pc + Xtip_gs;
991 double radius = fb==FRONT ? Cenx[0] - Xtip_pb/2 : Cenx[Nbrt] + Xtipt/2;
993 double Zhalfb = fb==FRONT ? (Bar_Z_cut-
Zmin)/2. : (
Zmax-
Zmin)/2.;
996 std::string
name = baseName +
"FrontBack::Absorber";
997 GeoBox *box =
new GeoBox(Xhalfb,Thce/2,Zhalfb);
998 GeoBox *box2 =
new GeoBox(Xhalfb,Thce/2,dz01);
999 const GeoLogVol *logVol =
new GeoLogVol(
name,box,Thin_abs);
1000 const GeoLogVol *logVol2 =
new GeoLogVol(
name,box2,Thick_abs);
1001 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
1002 GeoIntrusivePtr<GeoPhysVol> physVol2 =
new GeoPhysVol(logVol2);
1003 physVol->add(
new GeoTransform(GeoTrf::Translate3D(0.,0.,dz01-Zhalfb)));
1004 physVol->add(physVol2);
1005 GeoGenfun::GENFUNCTION Xcd =
radius*Cos(Gama);
1006 GeoGenfun::GENFUNCTION Ycd =
radius*Sin(Gama);
1007 GeoGenfun::GENFUNCTION Zcd = GeoGenfun::FixedConstant(
Zmin+Zhalfb);
1010 GeoXF::TRANSFUNCTION TX =
1011 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),Xcd) *
1012 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),Ycd) *
1013 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),Zcd) *
1014 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Gama);
1015 GeoSerialTransformer *st =
new GeoSerialTransformer(physVol, &TX, Nabsorber);
1016 m_ecamPhysicalPos->add(st);
1017 m_ecamPhysicalNeg->add(st);
1019 if (fb==FRONT) std::cout <<
" *** Front tip Absorber " << std::endl;
1020 else std::cout <<
" *** Back tip Absorber " << std::endl;
1021 std::cout <<
" Thin Abs Box " << Xhalfb <<
" " << Thce/2. <<
" " << Zhalfb << std::endl;
1022 std::cout <<
" Thick Abs Box " << Xhalfb <<
" " << Thce/2. <<
" " << dz01 << std::endl;
1023 std::cout <<
" Z position thick in thin " << dz01-Zhalfb << std::endl;
1024 std::cout <<
" Radial position " <<
radius << std::endl;
1025 std::cout <<
" Phi0 (Gaudi::Units::deg) " << Gama(0)/
Gaudi::Units::deg << std::endl;
1026 std::cout <<
" Z position in ECAM " <<
Zmin+Zhalfb << std::endl;
1033 double radiusg = Cenx[0]-Xtip_pb/2. - Xtips/2 ;
1034 double Zhalfbg = (Bar_Z_cut-
Zmin)/2. ;
1035 std::string
name = baseName +
"FrontBack::G10";
1037 std::cout <<
" *** Front tip G10 " << std::endl;
1038 std::cout <<
" Box parameters " << Xhalfbg <<
" " << Thce/2. <<
" " << Zhalfbg << std::endl;
1040 GeoBox *box =
new GeoBox(Xhalfbg,Thce/2,Zhalfbg);
1041 const GeoLogVol *logVol =
new GeoLogVol(
name,box,G10_bar);
1042 GeoPhysVol *physVol =
new GeoPhysVol(logVol);
1044 #ifdef BUILD_FRONT_STEEL
1047 name = baseName +
"FrontBack::Steel";
1048 double Tgfe = m_parameters->GetValue(
"LArEMBThinAbsIron");
1050 std::cout <<
" Iron Box parameters " << Xhalfbg
1051 <<
" " << Tgfe/4. <<
" " << Zhalfbg << std::endl;
1053 GeoBox *box2 =
new GeoBox(Xhalfbg,Tgfe/4.,Zhalfbg);
1054 const GeoLogVol *logVol2 =
new GeoLogVol(
name,box2,Iron);
1055 GeoIntrusivePtr<GeoPhysVol>physVol2 =
new GeoPhysVol(logVol2);
1057 std::cout <<
" Position Iron in G10 at y = +- " << 0.5*(+Thce-Tgfe/2.) << std::endl;
1059 physVol->add(
new GeoTransform(GeoTrf::Translate3D(0.,0.5*(-Thce+Tgfe/2.),0.)));
1060 physVol->add(physVol2);
1061 physVol->add(
new GeoTransform(GeoTrf::Translate3D(0.,0.5*(+Thce-Tgfe/2.),0.)));
1062 physVol->add(physVol2);
1063 #endif // build_front_steel
1066 GeoGenfun::GENFUNCTION Xcd = radiusg*Cos(Gama);
1067 GeoGenfun::GENFUNCTION Ycd = radiusg*Sin(Gama);
1068 GeoGenfun::GENFUNCTION Zcd = GeoGenfun::FixedConstant(Zhalfbg+
Zmin);
1069 GeoXF::TRANSFUNCTION TX =
1070 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),Xcd) *
1071 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),Ycd) *
1072 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),Zcd) *
1073 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Gama);
1074 GeoSerialTransformer *st =
new GeoSerialTransformer(physVol, &TX, Nabsorber);
1075 m_ecamPhysicalPos->add(st);
1076 m_ecamPhysicalNeg->add(st);
1078 std::cout <<
" Radial position G10 tip " << radiusg << std::endl;
1079 std::cout <<
" Phi0 (Gaudi::Units::deg)" << Gama(0)/
Gaudi::Units::deg << std::endl;
1080 std::cout <<
" Zposition in ECAM " <<
Zmin+Zhalfbg << std::endl;
1087 double Zhalfbe = fb==FRONT ? (Bar_Z_cut-
Zmin)/2. : (
Zmax -
Zmin)/2;
1088 double radiuse = fb==FRONT ? Cenx[0] - Xtips/2 : Cenx[Nbrt] + Xtipt/2;
1089 std::string
name = baseName +
"FrontBack::Electrode";
1090 GeoBox *box =
new GeoBox(Xhalfbe,Thel/2,Zhalfbe);
1091 const GeoLogVol *logVol =
new GeoLogVol(
name,box,Kapton_Cu);
1092 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
1093 GeoGenfun::GENFUNCTION Xcd = radiuse*Cos(Game);
1094 GeoGenfun::GENFUNCTION Ycd = radiuse*Sin(Game);
1095 GeoGenfun::GENFUNCTION Zcd = GeoGenfun::FixedConstant(
Zmin+Zhalfbe);
1097 GeoXF::TRANSFUNCTION TX =
1098 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),Xcd) *
1099 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),Ycd) *
1100 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),Zcd) *
1101 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Game);
1103 GeoSerialTransformer *st =
new GeoSerialTransformer(physVol, &TX, Nelectrode);
1104 m_ecamPhysicalPos->add(st);
1105 m_ecamPhysicalNeg->add(st);
1108 if (fb==FRONT) std::cout <<
" *** Front tip electrode " << std::endl;
1109 else std::cout <<
" *** Back tip electrode " << std::endl;
1110 std::cout <<
" Box " << Xhalfbe <<
" " << Thel/2. <<
" " << Zhalfbe << std::endl;
1111 std::cout <<
" Radial position " << radiuse << std::endl;
1112 std::cout <<
" Phi0 (Gaudi::Units::deg)" << Game(0)/
Gaudi::Units::deg << std::endl;
1113 std::cout <<
" Z position in ECAM " <<
Zmin+Zhalfbe << std::endl;
1131 for(
int irl=0; irl<Nbrt; irl++)
1134 std::cout <<
" ----- irl = " << irl << std::endl;
1138 if (checkParity==0) {
1139 iparit= (1-2*(irl%2));
1141 iparit=(2*(irl%2)-1);
1145 double Zcon1 = Zcp2[irl];
1146 double Zcon2 = Zcp2[irl+1];
1150 if (irl>0 && Rhoh[irl-1] < Bar_Rcmx) Zcon1=Zcp2h[irl-1];
1151 if (Rhoh[irl] < Bar_Rcmx) Zcon2=Zcp2h[irl];
1162 if (Rhocen[irl] < Bar_Rcmx && Rhocen[irl+1] > Bar_Rcmx) {
1164 std::cout <<
" Divide section in two pieces at r= " << Bar_Rcmx << std::endl;
1169 for (
int idivi=0;idivi<ndivi;idivi++) {
1183 frac=(Bar_Rcmx-Rhol[irl])/(Rhoh[irl]-Rhol[irl]);
1188 frac=(Rhoh[irl]-Bar_Rcmx)/(Rhoh[irl]-Rhol[irl]);
1191 std::cout <<
" Division " << idivi <<
" fraction of straight section " <<
frac << std::endl;
1195 tl1=tl1-safety_zlen;
1196 bl1=bl1-safety_zlen;
1202 double Dz = Thce/2.;
1205 GeoGenfun::GENFUNCTION x1a =
LArGeo::Fx(irl+0., Gama, Cenx, Ceny)
1208 GeoGenfun::GENFUNCTION x2a =
LArGeo::Fx(irl+1., Gama, Cenx, Ceny)
1211 GeoGenfun::GENFUNCTION y1a =
LArGeo::Fy(irl+0., Gama, Cenx, Ceny)
1214 GeoGenfun::GENFUNCTION y2a =
LArGeo::Fy(irl+1., Gama, Cenx, Ceny)
1217 GeoGenfun::GENFUNCTION
dx = x2a - x1a;
1218 GeoGenfun::GENFUNCTION
dy = y2a - y1a;
1222 GeoGenfun::GENFUNCTION Da = Sqrt (
dx*
dx +
dy*
dy );
1223 GeoGenfun::GENFUNCTION da = Sqrt ( (Da - 2.*Rint)*(Da + 2.*Rint) );
1226 GeoGenfun::GENFUNCTION cosalfa = (da*
dx -iparit*2.*Rint*
dy)/Da/Da;
1227 GeoGenfun::GENFUNCTION sinalfa = (da*
dy +iparit*2.*Rint*
dx)/Da/Da;
1228 GeoGenfun::GENFUNCTION newalpha =
LArGeo::ATan2(sinalfa,cosalfa);
1232 double Zx0 = (tl1+bl1)/2.;
1238 Xb1 = (Zcp1l[irl]-
Zmin)/2.;
1239 Xt1 = (Zcp1h[irl]-
Zmin)/2.;
1241 double xfrac=(Bar_Rcmx-Rhol[irl])/(Rhoh[irl]-Rhol[irl]);
1243 Xb1 = (Zcp1l[irl]-
Zmin)/2.;
1244 Xt1 = (xfrac*Zcp1h[irl]+(1.-xfrac)*Zcp1l[irl]-
Zmin)/2.;
1247 Xb1 = (xfrac*Zcp1h[irl]+(1.-xfrac)*Zcp1l[irl]-
Zmin)/2.;
1248 Xt1 = (Zcp1h[irl]-
Zmin)/2.;
1255 double Xt2 = tl1-Xt1;
1256 double Xb2 = bl1-Xb1;
1259 double Xtrans2 = Zx0 - (Xb2+Xt2)/2.;
1264 GeoGenfun::GENFUNCTION alpha = ATan(0.5*(bl1-tl1)/
h1);
1265 GeoGenfun::GENFUNCTION alpha_t = ATan(0.5*(Xb1-Xt1)/
h1);
1277 GeoGenfun::GENFUNCTION alpha_2 = ATan((2.*bl1-Xb2-(2.*tl1-Xt2))/(2.*
h1));
1283 GeoGenfun::GENFUNCTION alfrot = -
M_PI/2. - newalpha;
1285 GeoGenfun::GENFUNCTION Xcd = (x1a + x2a)/2. + (2.*idivi-1.)*(1.-
frac)*da/2.*cosalfa;
1286 GeoGenfun::GENFUNCTION Ycd = (y1a + y2a)/2. + (2.*idivi-1.)*(1.-
frac)*da/2.*sinalfa;
1287 GeoGenfun::GENFUNCTION Zcd = GeoGenfun::FixedConstant(
Zmin+(tl1+bl1)/2.+safety_zlen);
1289 GeoXF::TRANSFUNCTION TX =
1290 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),Xcd) *
1291 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),Ycd) *
1292 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),Zcd) *
1293 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),-alfrot)*
1302 std::string thinName = baseName +
"ThinAbs::Straight";
1303 std::string thickName = baseName +
"ThickAbs::Straight";
1306 std::cout <<
" *** thinAbs trap parameters " << std::endl;
1307 std::cout <<
" Thickness " <<
Dz << std::endl;
1308 std::cout <<
" h1 (rad leng) " <<
h1(
instance) << std::endl;
1309 std::cout <<
" tl1,bl1 " << tl1 <<
" " << bl1 << std::endl;
1310 std::cout <<
" alpha " << alpha(
instance) << std::endl;
1311 std::cout <<
" *** thickAbs trap parameters " << std::endl;
1312 std::cout <<
" Thickness " <<
Dz << std::endl;
1313 std::cout <<
" h1 (rad leng) " <<
h1(
instance) << std::endl;
1314 std::cout <<
" tl1,bl1 " << Xt1 <<
" " << Xb1 << std::endl;
1315 std::cout <<
" alpha " << alpha_t(
instance) << std::endl;
1317 std::cout <<
" Position absorber x,y,z " << Xcd(
instance) <<
" " << Ycd(
instance)
1318 <<
" " << Zcd(
instance) << std::endl;
1319 std::cout <<
" Rotation along Z " << -alfrot(
instance)/
deg << std::endl;
1320 std::cout <<
" Rotation along Y " << -90 << std::endl;
1325 GeoIntrusivePtr<GeoPhysVol> thinPhys =
nullptr;
1328 if (!doDetailedAbsorberStraight) {
1330 const GeoLogVol* thinLog =
new GeoLogVol(thinName,thinTrap,Thin_abs);
1331 thinPhys =
new GeoPhysVol(thinLog);
1335 const GeoLogVol* thickLog =
new GeoLogVol(thickName,thickTrap,Thick_abs);
1336 GeoIntrusivePtr<GeoPhysVol> thickPhys =
new GeoPhysVol(thickLog);
1338 thinPhys->add(
new GeoTransform(GeoTrf::TranslateX3D(Xtrans)));
1339 thinPhys->add(thickPhys);
1342 std::cout <<
" Position thick Abs in Thin at " << Xtrans << std::endl;
1346 if (doDetailedAbsorberStraight) {
1350 std::cout <<
" dz overall absorber volume " <<
Dz <<
" radial length " <<
h1(
instance) <<
1351 " lenghts " << tl1 <<
" " << bl1 <<
" Angle y axis " << alpha(
instance) << std::endl;
1353 const GeoLogVol* thinLog =
new GeoLogVol(thinName,thinTrap,myIron);
1354 thinPhys =
new GeoPhysVol(thinLog);
1358 double dz_glue =
Dz - 0.5*Thfe*Contract;
1360 std::cout <<
" dz glue volume in which lead will be put " << dz_glue << std::endl;
1362 GeoTrap* thickTrapGlue =
new GeoTrap(dz_glue,0.,0.,
h1(
instance),tl1,bl1,alpha(
instance),
1364 std::string thickGlueName = baseName +
"ThickAbsGlue::Straight";
1365 const GeoLogVol* thickTrapGlueLog =
new GeoLogVol(thickGlueName,thickTrapGlue, Glue);
1366 GeoIntrusivePtr<GeoPhysVol> thickTrapGluePhys =
new GeoPhysVol(thickTrapGlueLog);
1367 thinPhys->add(
new GeoTransform(GeoTrf::Translate3D(0.,0.,0.)));
1368 thinPhys->add(thickTrapGluePhys);
1371 double dz_lead_thick=0.5*Thpb*Contract;
1373 std::cout <<
" dz thick lead " << dz_lead_thick <<
" lenghts " << Xt1 <<
" " << Xb1 <<
" Angle y axis " << alpha_t(
instance) << std::endl;
1374 std::cout <<
" position in overall absorber at X trans " << Xtrans << std::endl;
1376 GeoTrap* thickTrapLead =
new GeoTrap(dz_lead_thick,0.,0.,
h1(
instance),Xt1,Xb1,alpha_t(
instance),
1378 std::string thickLeadName= baseName+
"ThickAbsLead::Straight";
1379 const GeoLogVol* thickTrapLeadLog =
new GeoLogVol(thickLeadName,thickTrapLead, myLead);
1380 GeoIntrusivePtr<GeoPhysVol> thickTrapLeadPhys =
new GeoPhysVol(thickTrapLeadLog);
1381 thickTrapGluePhys->add(
new GeoTransform(GeoTrf::TranslateX3D(Xtrans)));
1382 thickTrapGluePhys->add(thickTrapLeadPhys);
1384 double dz_lead_thin = 0.5*Thpb_thin*Contract;
1386 std::cout <<
" dz thin lead " << dz_lead_thin <<
" lenghts " << Xt2 <<
" " << Xb2 <<
" Angle y axis " << alpha_2(
instance) << std::endl;
1387 std::cout <<
" position in overall absorber at X trans " << Xtrans2 << std::endl;
1389 GeoTrap* thinTrapLead =
new GeoTrap(dz_lead_thin,0.,0.,
h1(
instance),Xt2,Xb2,alpha_2(
instance),
1391 std::string thinLeadName = baseName+
"ThinAbsLead::Straight";
1392 const GeoLogVol* thinTrapLeadLog =
new GeoLogVol(thinLeadName,thinTrapLead, myLead);
1393 GeoIntrusivePtr<GeoPhysVol> thinTrapLeadPhys =
new GeoPhysVol(thinTrapLeadLog);
1394 thickTrapGluePhys->add(
new GeoTransform(GeoTrf::TranslateX3D(Xtrans2)));
1395 thickTrapGluePhys->add(thinTrapLeadPhys);
1413 stacPhysical->add(
new GeoTransform(TX(
instance)));
1414 int icopytot=1000000*idivi+10000*irl+
instance;
1415 stacPhysical->add(
new GeoIdentifierTag(icopytot));
1416 stacPhysical->add(thinPhys);
1421 gStraightAbsorbers->
setHalfLength(irl,thinTrap->getDydzn());
1422 GeoSerialTransformer *st =
new GeoSerialTransformer(thinPhys, &TX, Nabsorber);
1423 stacPhysical->add(
new GeoSerialIdentifier(irl*10000+idivi*1000000));
1424 stacPhysical->add(st);
1432 double phi0_safety=0.;
1433 if (irl==0) phi0_safety=0.003;
1439 int irl2 = (irl==Nbrt-1) ? irl+1 : irl;
1441 for (
int jrl=irl1; jrl<=irl2; jrl++) {
1446 GeoGenfun::GENFUNCTION x0a =
LArGeo::Fx(iirl, Gama, Cenx, Ceny)
1449 GeoGenfun::GENFUNCTION y0a =
LArGeo::Fy(iirl, Gama, Cenx, Ceny)
1452 GeoGenfun::GENFUNCTION dx0 = x1a - x0a;
1453 GeoGenfun::GENFUNCTION dy0 = y1a - y0a;
1455 GeoGenfun::GENFUNCTION Da0 = Sqrt ( dx0*dx0 + dy0*dy0 );
1456 GeoGenfun::GENFUNCTION da0 = Sqrt ( (Da0 - 2.*Rint)*(Da0 + 2.*Rint) );
1458 GeoGenfun::GENFUNCTION cosalfa0 = (da0*dx0 +iparit*2.*Rint*dy0)/Da0/Da0;
1459 GeoGenfun::GENFUNCTION sinalfa0 = (da0*dy0 -iparit*2.*Rint*dx0)/Da0/Da0;
1460 GeoGenfun::GENFUNCTION alpha_prev =
LArGeo::ATan2(sinalfa0,cosalfa0);
1463 if (jrl>0 && jrl<Nbrt) {
1464 std::cout <<
"x1,y1,x0,y0 " << x1a(0) <<
" " << y1a(0) <<
" "
1465 << x0a(0) <<
" " << y0a(0) << std::endl;
1466 std::cout <<
" newalpha, alpha_prev " << newalpha(0) <<
" "
1467 << alpha_prev(0) << std::endl;
1470 GeoGenfun::Mod Mod2Pi(2*
M_PI);
1472 GeoGenfun::GENFUNCTION phi0_dfold_0 =
1473 GeoGenfun::FixedConstant(
M_PI/2.+phi0_safety);
1474 GeoGenfun::GENFUNCTION dphi_dfold_0 = Mod2Pi(newalpha-phi0_safety - Gama);
1475 GeoGenfun::GENFUNCTION phi0_dfold_1 = Mod2Pi(
M_PI/2.+ alpha_prev - Gama);
1476 GeoGenfun::GENFUNCTION dphi_dfold_1 = Mod2Pi(newalpha-alpha_prev);
1477 GeoGenfun::GENFUNCTION phi0_dfold_2 = Mod2Pi(
M_PI/2.+ newalpha - Gama);
1478 GeoGenfun::GENFUNCTION dphi_dfold_2 = Mod2Pi(- newalpha + Gama);
1480 GeoGenfun::GENFUNCTION phi0_ufold_0 =
1481 Mod2Pi(
M_PI/2.+newalpha-Gama);
1482 GeoGenfun::GENFUNCTION dphi_ufold_0 = Mod2Pi(-newalpha+Gama-phi0_safety);
1483 GeoGenfun::GENFUNCTION phi0_ufold_1 = Mod2Pi(
M_PI/2. + newalpha - Gama);
1484 GeoGenfun::GENFUNCTION dphi_ufold_1 = Mod2Pi(alpha_prev - newalpha);
1485 GeoGenfun::GENFUNCTION phi0_ufold_2 = GeoGenfun::FixedConstant(
M_PI/2.);
1486 GeoGenfun::GENFUNCTION dphi_ufold_2 = Mod2Pi(newalpha-Gama);
1488 const GeoGenfun::AbsFunction* phi0_fold=
nullptr;
1489 const GeoGenfun::AbsFunction* dphi_fold=
nullptr;
1490 const GeoXF::Function* TXfold=
nullptr;
1492 std::string thinName;
1493 std::string thickName;
1494 if (jrl%2==checkParity) {
1495 thinName = baseName +
"ThinAbs::CornerDownFold";
1496 thickName = baseName +
"ThickAbs::CornerDownFold";
1498 thinName = baseName +
"ThinAbs::CornerUpFold";
1499 thickName = baseName +
"ThickAbs::CornerUpFold";
1503 double Rcmin=Rint-Thce/2.;
1504 double Rcmax=Rint+Thce/2.;
1507 double ddz0 = dz0-safety_zlen;
1508 double ddz01 = dz01-safety_zlen;
1510 ddz0=dza-safety_zlen;
1511 ddz01=dza1-safety_zlen;
1515 GeoGenfun::GENFUNCTION zpos = GeoGenfun::FixedConstant(
Zmin+dz0);
1517 if (jrl%2==checkParity) phirot = -
M_PI;
1518 GeoXF::TRANSFUNCTION TXfold1=
1519 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),x1a) *
1520 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),y1a) *
1521 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),zpos) *
1522 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Gama+phirot);
1523 GeoXF::TRANSFUNCTION TXfold2 =
1524 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),x2a) *
1525 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),y2a) *
1526 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),zpos) *
1527 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Gama+phirot);
1530 if (jrl==0 && checkParity==0) {
1531 phi0_fold = &(phi0_dfold_0);
1532 dphi_fold = &(dphi_dfold_0);
1533 TXfold = &(TXfold1);
1536 if (jrl%2==checkParity && jrl >0 && jrl<Nbrt) {
1537 phi0_fold = &(phi0_dfold_1);
1538 dphi_fold = &(dphi_dfold_1);
1539 TXfold = &(TXfold1);
1542 if (jrl%2 ==checkParity && jrl == Nbrt) {
1543 phi0_fold = &(phi0_dfold_2);
1544 dphi_fold = &(dphi_dfold_2);
1545 TXfold = &(TXfold2);
1548 if (jrl==0 && checkParity==1) {
1549 phi0_fold = &(phi0_ufold_0);
1550 dphi_fold = &(dphi_ufold_0);
1551 TXfold = &(TXfold1);
1554 if (jrl%2 ==(1-checkParity) && jrl>0 && jrl < Nbrt) {
1555 phi0_fold = &(phi0_ufold_1);
1556 dphi_fold = &(dphi_ufold_1);
1557 TXfold = &(TXfold1);
1560 if (jrl%2==(1-checkParity) && jrl==Nbrt) {
1561 phi0_fold = &(phi0_ufold_2);
1562 dphi_fold = &(dphi_ufold_2);
1563 TXfold = &(TXfold2);
1566 if (!phi0_fold || !dphi_fold || !TXfold) {
1567 log << MSG::INFO <<
" LArGeoBarrel::BarrelConstruction fold not defined..." <<
endmsg;
1573 GeoIntrusivePtr<GeoPhysVol> thinPhys =
nullptr;
1575 if (!doDetailedAbsorberFold) {
1576 GeoTubs* thinTubs =
new GeoTubs(Rcmin,Rcmax,ddz0,
1578 GeoTubs* thickTubs =
new GeoTubs(Rcmin,Rcmax,ddz01,
1581 std::cout <<
" Thin Abs Corner (tubs) " << std::endl;
1582 std::cout <<
" Rmin,Rmax,dZ,phi0,Dphi[rad] " << Rcmin <<
" " << Rcmax <<
" "
1583 << ddz0 <<
" " << (*phi0_fold)(
instance) <<
" "
1584 <<(*dphi_fold)(
instance) << std::endl;
1585 std::cout <<
" Thick Abs Corner (tubs) " << std::endl;
1586 std::cout <<
" Rmin,Rmax,dZ,phi0,Dphi[rad] " << Rcmin <<
" " << Rcmax <<
" "
1587 << ddz01 <<
" " << (*phi0_fold)(
instance) <<
" "
1588 << (*dphi_fold)(
instance) << std::endl;
1590 const GeoLogVol* thinLog =
new GeoLogVol(thinName,thinTubs,Thin_abs);
1591 const GeoLogVol* thickLog =
new GeoLogVol(thickName,thickTubs,Thick_abs);
1593 thinPhys =
new GeoPhysVol(thinLog);
1594 GeoIntrusivePtr<GeoPhysVol> thickPhys =
new GeoPhysVol(thickLog);
1596 thinPhys->add(
new GeoTransform(GeoTrf::TranslateZ3D(ddz01-ddz0)));
1597 thinPhys->add(thickPhys);
1599 std::cout <<
" Position Thick fold in Thin Z = " << ddz01-ddz0 << std::endl;
1604 if (doDetailedAbsorberFold) {
1607 GeoTubs* thinTubs =
new GeoTubs(Rcmin,Rcmax,ddz0,
1609 const GeoLogVol* thinLog =
new GeoLogVol(thinName,thinTubs,myIron);
1610 thinPhys =
new GeoPhysVol(thinLog);
1613 std::cout <<
" overall fold volume rmin,rmax,dz " << Rcmin <<
" " << Rcmax <<
" " << ddz0 << std::endl;
1617 GeoTubs* glueTubs =
new GeoTubs(Rcmin+0.5*Thfe*Contract,Rcmax-0.5*Thfe*Contract,ddz0,
1619 std::string foldGlueName = baseName+
"Glue::Fold";
1620 const GeoLogVol* glueTubsLog =
new GeoLogVol(foldGlueName,glueTubs,Glue);
1621 GeoIntrusivePtr<GeoPhysVol> glueTubsPhys =
new GeoPhysVol(glueTubsLog);
1622 thinPhys->add(
new GeoTransform(GeoTrf::TranslateZ3D(0.)));
1623 thinPhys->add(glueTubsPhys);
1625 std::cout <<
" glue fold volume " << Rcmin+0.5*Thfe*Contract <<
" " << Rcmax-0.5*Thfe*Contract <<
" " << ddz0 << std::endl;
1631 GeoTubs* thickLeadTubs =
new GeoTubs(Rint-0.5*Thpb*Contract,Rint+0.5*Thpb*Contract,ddz01,
1633 std::string foldThickLeadName = baseName+
"ThickLead::Fold";
1634 const GeoLogVol* thickLeadLog =
new GeoLogVol(foldThickLeadName,thickLeadTubs,myLead);
1635 GeoIntrusivePtr<GeoPhysVol> thickLeadPhys =
new GeoPhysVol(thickLeadLog);
1636 glueTubsPhys->add(
new GeoTransform(GeoTrf::TranslateZ3D(ddz01-ddz0)));
1637 glueTubsPhys->add(thickLeadPhys);
1639 std::cout <<
" thick lead volume " << Rint-Thpb*Contract <<
" " << Rint+Thpb*Contract <<
" " << ddz01 << std::endl;
1640 std::cout <<
" put in glue fold volume with z translation " << ddz01-ddz0 << std::endl;
1644 GeoTubs* thinLeadTubs =
new GeoTubs(Rint-0.5*Thpb_thin*Contract,Rint+0.5*Thpb_thin*Contract,ddz0-ddz01,
1646 std::string foldThinLeadName = baseName+
"ThinLead::Fold";
1647 const GeoLogVol* thinLeadLog =
new GeoLogVol(foldThinLeadName,thinLeadTubs,myLead);
1648 GeoIntrusivePtr<GeoPhysVol> thinLeadPhys =
new GeoPhysVol(thinLeadLog);
1649 glueTubsPhys->add(
new GeoTransform(GeoTrf::TranslateZ3D(ddz01)));
1650 glueTubsPhys->add(thinLeadPhys);
1653 std::cout <<
" thin lead volume " << Rint-Thpb_thin*Contract <<
" " << Rint+Thpb_thin*Contract <<
" " << ddz01 << std::endl;
1654 std::cout <<
" put in glue fold volume with z translation " << ddz01 << std::endl;
1663 std::cout <<
" Position absorber fold x,y,z,rotation/z " << x1a(
instance) <<
1667 std::cout <<
" Position absorber fold x,y,z,rotation/z " << x2a(
instance) <<
1673 stacPhysical->add(
new GeoTransform((*TXfold)(
instance)));
1675 stacPhysical->add(
new GeoIdentifierTag(icopytot));
1676 stacPhysical->add(thinPhys);
1679 GeoSerialTransformer *st =
new GeoSerialTransformer(thinPhys, TXfold, Nabsorber);
1680 stacPhysical->add(
new GeoSerialIdentifier(jrl*10000));
1681 stacPhysical->add(st);
1695 double Dze = Thel/2.;
1698 GeoGenfun::GENFUNCTION x1e =
LArGeo::Fx(irl+0., Game, Cenx, Ceny)
1701 GeoGenfun::GENFUNCTION x2e =
LArGeo::Fx(irl+1., Game, Cenx, Ceny)
1704 GeoGenfun::GENFUNCTION y1e =
LArGeo::Fy(irl+0., Game, Cenx, Ceny)
1707 GeoGenfun::GENFUNCTION y2e =
LArGeo::Fy(irl+1., Game, Cenx, Ceny)
1710 GeoGenfun::GENFUNCTION dxe = x2e - x1e;
1711 GeoGenfun::GENFUNCTION dye = y2e - y1e;
1713 GeoGenfun::GENFUNCTION De = Sqrt ( dxe*dxe + dye*dye );
1714 GeoGenfun::GENFUNCTION de = Sqrt ( (De - 2.*Rint)*(De + 2.*Rint) );
1717 GeoGenfun::GENFUNCTION cosalfae = (de*dxe -iparit*2.*Rint*dye)/De/De;
1718 GeoGenfun::GENFUNCTION sinalfae = (de*dye +iparit*2.*Rint*dxe)/De/De;
1719 GeoGenfun::GENFUNCTION newalphe =
LArGeo::ATan2(sinalfae,cosalfae);
1725 GeoGenfun::GENFUNCTION alfrote = -
M_PI/2. - newalphe;
1727 GeoGenfun::GENFUNCTION Xcde = (x1e + x2e)/2.+ (2.*idivi-1.)*(1.-
frac)*de/2.*cosalfae;
1728 GeoGenfun::GENFUNCTION Ycde = (y1e + y2e)/2.+ (2.*idivi-1.)*(1.-
frac)*de/2.*sinalfae;
1729 GeoGenfun::GENFUNCTION Zcde = GeoGenfun::FixedConstant(
Zmin+(tl1+bl1)/2.+safety_zlen);
1733 GeoGenfun::GENFUNCTION alpha_e = ATan(0.5*(bl1-tl1)/h1e);
1736 GeoXF::TRANSFUNCTION TXE =
1737 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),Xcde) *
1738 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),Ycde) *
1739 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),Zcde) *
1740 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),-alfrote)*
1746 std::string
name = baseName +
"Electrode::Straight";
1748 std::cout <<
" *** Electrode trap parameters " << std::endl;
1749 std::cout <<
" Thickness " << Dze << std::endl;
1750 std::cout <<
" h1 (rad leng) " << h1e(
instance) << std::endl;
1751 std::cout <<
" tl1,tb1 " << tl1 <<
" " << bl1 << std::endl;
1752 std::cout <<
" alpha " << alpha_e(
instance) << std::endl;
1753 std::cout <<
" Position electrode x,y,z " << Xcde(
instance) <<
" " << Ycde(
instance)
1754 <<
" " << Zcde(
instance) << std::endl;
1755 std::cout <<
" Rotation along Z " << -alfrote(
instance)/
deg << std::endl;
1756 std::cout <<
" Rotation along Y " << -90 << std::endl;
1760 GeoTrap* trap =
new GeoTrap(Dze,0.,0.,h1e(
instance),tl1,bl1,alpha_e(
instance),
1762 const GeoLogVol* logVol =
new GeoLogVol(
name,trap,Kapton_Cu);
1763 GeoIntrusivePtr<GeoPhysVol> physVol =
new GeoPhysVol(logVol);
1777 stacPhysical->add(
new GeoTransform(TXE(
instance)));
1778 int icopytot=1000000*idivi+10000*irl+
instance;
1779 stacPhysical->add(
new GeoIdentifierTag(icopytot));
1780 stacPhysical->add(physVol);
1786 GeoSerialTransformer *st =
new GeoSerialTransformer(physVol, &TXE, Nelectrode);
1787 stacPhysical->add(
new GeoSerialIdentifier(irl*10000+idivi*1000000));
1788 stacPhysical->add(st);
1797 double phi0_safety=0.;
1798 if (irl==0) phi0_safety=0.003;
1804 int irl2 = (irl==Nbrt-1) ? irl+1 : irl;
1806 for (
int jrl=irl1; jrl<=irl2; jrl++) {
1811 GeoGenfun::GENFUNCTION x0e =
LArGeo::Fx(iirl, Game, Cenx, Ceny)
1814 GeoGenfun::GENFUNCTION y0e =
LArGeo::Fy(iirl, Game, Cenx, Ceny)
1817 GeoGenfun::GENFUNCTION dx0 = x1e - x0e;
1818 GeoGenfun::GENFUNCTION dy0 = y1e - y0e;
1820 GeoGenfun::GENFUNCTION Da0 = Sqrt ( dx0*dx0 + dy0*dy0 );
1821 GeoGenfun::GENFUNCTION da0 = Sqrt ( (Da0 - 2.*Rint)*(Da0 + 2.*Rint) );
1823 GeoGenfun::GENFUNCTION cosalfa0 = (da0*dx0 +iparit*2.*Rint*dy0)/Da0/Da0;
1824 GeoGenfun::GENFUNCTION sinalfa0 = (da0*dy0 -iparit*2.*Rint*dx0)/Da0/Da0;
1825 GeoGenfun::GENFUNCTION alphe_prev =
LArGeo::ATan2(sinalfa0,cosalfa0);
1828 if (jrl>0 && jrl<Nbrt) {
1829 std::cout <<
" newalphae, alphae_prev " << newalphe(0) <<
" "
1830 << alphe_prev(0) << std::endl;
1835 GeoGenfun::Mod Mod2Pi(2*
M_PI);
1836 GeoGenfun::GENFUNCTION phi0_dfold_0 =
1837 GeoGenfun::FixedConstant(
M_PI/2.+phi0_safety);
1838 GeoGenfun::GENFUNCTION dphi_dfold_0 = Mod2Pi(newalphe-phi0_safety-Game);
1839 GeoGenfun::GENFUNCTION phi0_dfold_1 = Mod2Pi(
M_PI/2.+ alphe_prev - Game);
1840 GeoGenfun::GENFUNCTION dphi_dfold_1 = Mod2Pi(newalphe-alphe_prev);
1841 GeoGenfun::GENFUNCTION phi0_dfold_2 = Mod2Pi(
M_PI/2.+ newalphe - Game);
1842 GeoGenfun::GENFUNCTION dphi_dfold_2 = Mod2Pi(- newalphe + Game);
1844 GeoGenfun::GENFUNCTION phi0_ufold_0 =
1845 Mod2Pi(
M_PI/2.+newalphe-Game);
1846 GeoGenfun::GENFUNCTION dphi_ufold_0 = Mod2Pi(-newalphe+Game-phi0_safety);
1847 GeoGenfun::GENFUNCTION phi0_ufold_1 = Mod2Pi(
M_PI/2. + newalphe - Game);
1848 GeoGenfun::GENFUNCTION dphi_ufold_1 = Mod2Pi(alphe_prev - newalphe);
1849 GeoGenfun::GENFUNCTION phi0_ufold_2 = GeoGenfun::FixedConstant(
M_PI/2.);
1850 GeoGenfun::GENFUNCTION dphi_ufold_2 = Mod2Pi(newalphe - Game);
1852 const GeoGenfun::AbsFunction* phi0_fold=
nullptr;
1853 const GeoGenfun::AbsFunction* dphi_fold=
nullptr;
1854 const GeoXF::Function* TXfold=
nullptr;
1857 if (jrl%2==checkParity) {
1858 eName = baseName +
"Electrode::CornerDownFold";
1860 eName = baseName +
"Electrode::CornerUpFold";
1864 double Rcmine=Rint-Thel/2.;
1865 double Rcmaxe=Rint+Thel/2.;
1868 double ddz0 = dz0-safety_zlen;
1870 ddz0 = dza - safety_zlen;
1873 GeoGenfun::GENFUNCTION zpos = GeoGenfun::FixedConstant(
Zmin+dz0);
1875 if (jrl%2==checkParity) phirot = -
M_PI;
1876 GeoXF::TRANSFUNCTION TXfold1=
1877 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),x1e) *
1878 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),y1e) *
1879 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),zpos) *
1880 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Game+phirot);
1881 GeoXF::TRANSFUNCTION TXfold2 =
1882 GeoXF::Pow(GeoTrf::TranslateX3D(1.0),x2e) *
1883 GeoXF::Pow(GeoTrf::TranslateY3D(1.0),y2e) *
1884 GeoXF::Pow(GeoTrf::TranslateZ3D(1.0),zpos) *
1885 GeoXF::Pow(GeoTrf::RotateZ3D(1.0),Game+phirot);
1888 if (jrl==0 && checkParity==0) {
1889 phi0_fold = &(phi0_dfold_0);
1890 dphi_fold = &(dphi_dfold_0);
1891 TXfold = &(TXfold1);
1894 if (jrl%2==checkParity && jrl >0 && jrl<Nbrt) {
1895 phi0_fold = &(phi0_dfold_1);
1896 dphi_fold = &(dphi_dfold_1);
1897 TXfold = &(TXfold1);
1900 if (jrl%2 ==checkParity && jrl == Nbrt) {
1901 phi0_fold = &(phi0_dfold_2);
1902 dphi_fold = &(dphi_dfold_2);
1903 TXfold = &(TXfold2);
1906 if (jrl==0 && checkParity==1 ) {
1907 phi0_fold = &(phi0_ufold_0);
1908 dphi_fold = &(dphi_ufold_0);
1909 TXfold = &(TXfold1);
1913 if (jrl%2 ==(1-checkParity) && jrl>0 && jrl < Nbrt) {
1914 phi0_fold = &(phi0_ufold_1);
1915 dphi_fold = &(dphi_ufold_1);
1916 TXfold = &(TXfold1);
1919 if (jrl%2==(1-checkParity) && jrl==Nbrt) {
1920 phi0_fold = &(phi0_ufold_2);
1921 dphi_fold = &(dphi_ufold_2);
1922 TXfold = &(TXfold2);
1925 if (!phi0_fold || !dphi_fold || !TXfold) {
1926 log << MSG::INFO <<
" LArGeoBarrel::BarrelConstruction fold not defined..." <<
endmsg;
1932 GeoTubs* foldeTubs =
new GeoTubs(Rcmine,Rcmaxe,ddz0,
1935 std::cout <<
" Electrode Corner (tubs) " << std::endl;
1936 std::cout <<
" Rmin,Rmax,dZ,phi0,Dphi[rad] " << Rcmine <<
" " << Rcmaxe <<
" "
1937 << ddz0 <<
" " << (*phi0_fold)(
instance) <<
" "
1938 <<(*dphi_fold)(
instance) << std::endl;
1940 const GeoLogVol* foldeLog =
new GeoLogVol(eName,foldeTubs,Kapton_Cu);
1942 GeoIntrusivePtr<GeoPhysVol> foldePhys =
new GeoPhysVol(foldeLog);
1946 std::cout <<
" Position electrode fold x,y,z,rotation/z " << x1e(
instance) <<
1950 std::cout <<
" Position electrode fold x,y,z,rotation/z " << x2e(
instance) <<
1956 stacPhysical->add(
new GeoTransform((*TXfold)(
instance)));
1958 stacPhysical->add(
new GeoIdentifierTag(icopytot));
1959 stacPhysical->add(foldePhys);
1962 GeoSerialTransformer *st =
new GeoSerialTransformer(foldePhys, TXfold, Nelectrode);
1963 stacPhysical->add(
new GeoSerialIdentifier(jrl*10000));
1964 stacPhysical->add(st);
1977 if(!
status.isSuccess())
throw std::runtime_error (
"Cannot store STRAIGHTELECTRODES structure");
1981 if(!
status.isSuccess())
throw std::runtime_error (
"Cannot store STRAIGHTABSORBERS structure");
1985 #endif // BUILD_ACCORDION_PLATES
1987 log <<
MSG::DEBUG <<
"++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
1988 log <<
"+ +" << std::endl;
1989 log <<
"+ END of Barrel EM GeoModel definition +" << std::endl;
1990 log <<
"+ +" << std::endl;
1991 log <<
"++++++++++++++++++++++++++++++++++++++++++++++++++++" <<
endmsg;
1998 std::cout <<
"******** Print out of LArBarrel Accordion parameters " << std::endl;
1999 std::cout <<
"Absorber parameters: " <<std::endl;
2000 std::cout <<
"Glue for ThinAbs " <<
2001 m_parameters->GetValue(
"LArEMBThinAbsGlue") << std::endl;
2002 std::cout <<
"Iron for ThinAbs " <<
2003 m_parameters->GetValue(
"LArEMBThinAbsIron") << std::endl;
2004 std::cout <<
"Lead for ThinAbs " <<
2005 m_parameters->GetValue(
"LArEMBThinAbsLead") << std::endl;
2006 std::cout <<
"Glue for ThickAbs " <<
2007 m_parameters->GetValue(
"LArEMBThickAbsGlue") << std::endl;
2008 std::cout <<
"Iron for ThickAbs " <<
2009 m_parameters->GetValue(
"LArEMBThickAbsIron") << std::endl;
2010 std::cout <<
"Lead for ThickAbs "
2011 << m_parameters->GetValue(
"LArEMBThickAbsLead") << std::endl;
2013 std::cout <<
"Electrode parameters: " << std::endl;
2014 std::cout <<
"Copper thickness " <<
2015 m_parameters->GetValue(
"LArEMBThickElecCopper") << std::endl;
2016 std::cout <<
"Kapton thickness " <<
2017 m_parameters->GetValue(
"LArEMBThickElecKapton") << std::endl;
2019 std::cout <<
"Accordion geometry description " << std::endl;
2020 std::cout <<
"phi first absorber " <<
2021 m_parameters->GetValue(
"LArEMBAbsPhiFirst") << std::endl;
2022 std::cout <<
"zmin of absorber " <<
2023 m_parameters->GetValue(
"LArEMBMotherZmin") << std::endl;
2024 std::cout <<
"zmax of absorber " <<
2025 m_parameters->GetValue(
"LArEMBMotherZmax") << std::endl;
2026 std::cout <<
"Max eta " <<
2027 m_parameters->GetValue(
"LArEMBMaxEtaAcceptance") << std::endl;
2028 std::cout <<
"Nominal eta lead change " <<
2029 m_parameters->GetValue(
"LArEMBThickEtaAcceptance") << std::endl;
2030 std::cout <<
"Number of abs in phi " <<
2031 m_parameters->GetValue(
"LArEMBnoOFPhysPhiCell") << std::endl;
2032 std::cout <<
"Dphi between absorbers " <<
2033 m_parameters->GetValue(
"LArEMBPhiGapAperture") << std::endl;
2034 std::cout <<
"Rmin Mother " <<
2035 m_parameters->GetValue(
"LArEMBMotherRmin") << std::endl;
2036 std::cout <<
"Rmax Mother " <<
2037 m_parameters->GetValue(
"LArEMBMotherRmax") << std::endl;
2038 std::cout <<
"Rmin at zmax " <<
2039 m_parameters->GetValue(
"LArEMBRminHighZ") << std::endl;
2040 std::cout <<
"Phimax barrel " <<
2041 m_parameters->GetValue(
"LArEMBphiMaxBarrel") << std::endl;
2042 std::cout <<
"Number of zigs " <<
2043 m_parameters->GetValue(
"LArEMBnoOFAccZigs") << std::endl;
2044 std::cout <<
"Fold Gaudi::Units::rad of curvature " <<
2045 m_parameters->GetValue(
"LArEMBNeutFiberRadius") << std::endl;
2046 for (
int i=0;
i<15;
i++) {
2047 std::cout <<
"Fold " <<
i <<
" radius " <<
2048 m_parameters->GetValue(
"LArEMBRadiusAtCurvature",
i) << std::endl;
2050 for (
int i=0;
i<15;
i++) {
2051 std::cout <<
"Fold " <<
i <<
" phi0 " <<
2052 m_parameters->GetValue(
"LArEMBPhiAtCurvature",
i) << std::endl;
2054 for (
int i=0;
i<15;
i++) {
2055 std::cout <<
"Fold " <<
i <<
" eta trans " <<
2056 m_parameters->GetValue(
"LArEMBEtaTrans",
i) << std::endl;
2058 for (
int i=0;
i<14;
i++) {
2059 std::cout <<
"zig " <<
i <<
" angle " <<
2060 m_parameters->GetValue(
"LArEMBDeltaZigAngle",
i) << std::endl;
2063 std::cout <<
"Fiducial active region definition " << std::endl;
2064 std::cout <<
"Zmin " <<
2065 m_parameters->GetValue(
"LArEMBfiducialMothZmin") << std::endl;
2066 std::cout <<
"Zmax " <<
2067 m_parameters->GetValue(
"LArEMBfiducialMothZmax") << std::endl;
2068 std::cout <<
"Rmin " <<
2069 m_parameters->GetValue(
"LArEMBRadiusInnerAccordion") << std::endl;
2070 std::cout <<
"Rmax " <<
2071 m_parameters->GetValue(
"LArEMBFiducialRmax") << std::endl;
2072 std::cout <<
"Inactive DR S1 to S2 " <<
2073 m_parameters->GetValue(
"LArEMBDeltaRS12") << std::endl;
2075 std::cout <<
"Description of matter from PS to strip " << std::endl;
2076 std::cout <<
"DeltaR for LAR+Cable+MotherBoard volume " <<
2077 m_parameters->GetValue(
"LArEMBInnerElectronics") << std::endl;
2078 std::cout <<
"Number of cable volumes " <<
2079 m_parameters->GetValue(
"LArEMBnoOFcableBundle") << std::endl;
2080 std::cout <<
"Cable Cu fraction " <<
2081 m_parameters->GetValue(
"LArEMBmasspercentCu") << std::endl;
2082 std::cout <<
"Cable Kapton fraction " <<
2083 m_parameters->GetValue(
"LArEMBmasspercentKap") << std::endl;
2084 std::cout <<
"Cable thickness@eta=0 " <<
2085 m_parameters->GetValue(
"LArEMBCablethickat0") << std::endl;
2086 std::cout <<
"Cable thick increase/eta " <<
2087 m_parameters->GetValue(
"LArEMBthickincrfac") << std::endl;
2088 std::cout <<
"Cable width " <<
2089 m_parameters->GetValue(
"LArEMBCableEtaheight") << std::endl;
2090 std::cout <<
"DR of cable from PS " <<
2091 m_parameters->GetValue(
"LArEMBCablclearfrPS") << std::endl;
2092 std::cout <<
"Number of motherboard " <<
2093 m_parameters->GetValue(
"LArEMBnoOFmothboard") << std::endl;
2094 std::cout <<
"MotherBoard Cu thickness " <<
2095 m_parameters->GetValue(
"LArEMBCuThickness") << std::endl;
2096 std::cout <<
"MotherBoard G10 thickness " <<
2097 m_parameters->GetValue(
"LArEMBG10Thickness") << std::endl;
2098 std::cout <<
"MotherBoard thickness " <<
2099 m_parameters->GetValue(
"LArEMBMoBoTchickness") << std::endl;
2100 std::cout <<
"MotherBoard width " <<
2101 m_parameters->GetValue(
"LArEMBMoBoHeight") << std::endl;
2102 std::cout <<
"MotherBoard DR from PS " <<
2103 m_parameters->GetValue(
"LArEMBMoBoclearfrPS") << std::endl;
2104 std::cout <<
"G10 inner ring thickness " <<
2105 m_parameters->GetValue(
"LArEMBG10SupportBarsIn") << std::endl;
2106 std::cout <<
"Dz G10 inner ring " <<
2107 m_parameters->GetValue(
"LArEMBG10FrontDeltaZ") << std::endl;
2108 std::cout <<
"G10 front tip " <<
2109 m_parameters->GetValue(
"LArEMBG10TipThickFront") << std::endl;
2110 std::cout <<
"Absorber front tip " <<
2111 m_parameters->GetValue(
"LArEMBLeadTipThickFront") << std::endl;
2112 std::cout <<
"Description of matter after accordion " << std::endl;
2113 std::cout <<
"G10 outer ring thickness " <<
2114 m_parameters->GetValue(
"LArEMBG10SupportBarsOut") << std::endl,
2115 std::cout <<
"Absorber outer tip " <<
2116 m_parameters->GetValue(
"LArEMBLeadTipThickEnd") << std::endl;
2117 std::cout <<
"G10 outer tip " <<
2118 m_parameters->GetValue(
"LArEMBG10TipThickEnd") << std::endl;
2119 std::cout <<
"total outer tip DeltaR " <<
2120 m_parameters->GetValue(
"LArEMBLArGapTail") << std::endl;