diff --git a/CalPatRec/src/CalHelixFinder_module.cc b/CalPatRec/src/CalHelixFinder_module.cc index 697e5f5de6..b7ec9a973b 100644 --- a/CalPatRec/src/CalHelixFinder_module.cc +++ b/CalPatRec/src/CalHelixFinder_module.cc @@ -67,8 +67,8 @@ namespace mu2e { std::vector helvals = config().Helicities(); for(auto hv : helvals) { - Helicity hel(hv); - _hels.push_back(hel); + Helicity hel(hv); + _hels.push_back(hel); } if (_doSingleOutput){ produces(); @@ -79,9 +79,9 @@ namespace mu2e { } _tpart = (TrkParticle::type)_fitparticle; -//----------------------------------------------------------------------------- -// provide for interactive disanostics -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // provide for interactive disanostics + //----------------------------------------------------------------------------- _helTraj = 0; _data.shLabel = _shLabel; @@ -92,20 +92,20 @@ namespace mu2e { else _hmanager = std::make_unique(); } -//----------------------------------------------------------------------------- -// destructor -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // destructor + //----------------------------------------------------------------------------- CalHelixFinder::~CalHelixFinder() { if (_helTraj) delete _helTraj; } -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- void CalHelixFinder::beginJob(){ art::ServiceHandle tfs; _hmanager->bookHistograms(tfs); } -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- void CalHelixFinder::beginRun(art::Run& ) { mu2e::GeomHandle bfmgr; mu2e::GeomHandle det; @@ -168,9 +168,9 @@ namespace mu2e { } -//----------------------------------------------------------------------------- -// find the input data objects -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // find the input data objects + //----------------------------------------------------------------------------- bool CalHelixFinder::findData(const art::Event& evt) { if (evt.getByLabel(_shLabel, _strawhitsH)) { @@ -179,7 +179,7 @@ namespace mu2e { else { _chcol = nullptr; printf(" >>> ERROR in CalHelixFinder::findData: StrawHitCollection with label=%s not found.\n", - _shLabel.data()); + _shLabel.data()); } // art::Handle shposH; @@ -198,22 +198,22 @@ namespace mu2e { else { _timeclcol = nullptr; printf(" >>> ERROR in CalHelixFinder::findData: TimeClusterCollection with label=%s not found.\n", - _timeclLabel.data()); + _timeclLabel.data()); } -//----------------------------------------------------------------------------- -// done -//----------------------------------------------------------------------------- - return (_chcol != nullptr) //&& (_shpcol != nullptr) - && (_timeclcol != nullptr); + //----------------------------------------------------------------------------- + // done + //----------------------------------------------------------------------------- + return (_chcol != nullptr) //&& (_shpcol != nullptr) + && (_timeclcol != nullptr); } -//----------------------------------------------------------------------------- -// event entry point -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // event entry point + //----------------------------------------------------------------------------- void CalHelixFinder::produce(art::Event& event ) { const char* oname = "CalHelixFinder::filter"; - // diagnostic info + // diagnostic info _data.event = &event; _iev = event.id().event(); int nGoodTClusterHits(0); @@ -233,16 +233,16 @@ namespace mu2e { _data.nseeds [counter] = 0; } // unique_ptr outseeds(new HelixSeedCollection); -//----------------------------------------------------------------------------- -// find the data -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // find the data + //----------------------------------------------------------------------------- if (!findData(event)) { printf("%s ERROR: No straw hits found, RETURN\n", oname); - goto END; + goto END; } -//----------------------------------------------------------------------------- -// loop over found time peaks - for us, - "eligible" calorimeter clusters -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // loop over found time peaks - for us, - "eligible" calorimeter clusters + //----------------------------------------------------------------------------- _hfResult._tpart = _tpart; _hfResult._fdir = _fdir; _hfResult._chcol = _chcol; @@ -258,28 +258,28 @@ namespace mu2e { std::vector helix_seed_vec; -//----------------------------------------------------------------------------- -// create track definitions for the helix fit from this initial information -// track fitting objects for this peak -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // create track definitions for the helix fit from this initial information + // track fitting objects for this peak + //----------------------------------------------------------------------------- _hfResult.clearTempVariables();//clearTimeClusterInfo(); _hfResult._timeCluster = tc; _hfResult._timeClusterPtr = art::Ptr(_timeclcolH,ipeak); -//----------------------------------------------------------------------------- -// fill the face-order hits collector -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // fill the face-order hits collector + //----------------------------------------------------------------------------- _hfinder.fillFaceOrderedHits(_hfResult); -//----------------------------------------------------------------------------- -// Step 1: now loop over the two possible helicities. -// Find initial helical approximation of a track for both hypothesis -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // Step 1: now loop over the two possible helicities. + // Find initial helical approximation of a track for both hypothesis + //----------------------------------------------------------------------------- std::vector nHitsRatio_vec; for (size_t i=0; i<_hels.size(); ++i){ -//----------------------------------------------------------------------------- -// create track definitions for the helix fit from this initial information -// track fitting objects for this peak -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // create track definitions for the helix fit from this initial information + // track fitting objects for this peak + //----------------------------------------------------------------------------- CalHelixFinderData tmpResult(_hfResult); tmpResult.clearHelixInfo(); @@ -290,7 +290,7 @@ namespace mu2e { if (!rc) continue; if(!tmpResult.helix()) { std::cout << "[CalHelixFinder::" << __func__ << "] " << event.id() - << " Helix found but nullptr returned!!\n"; + << " Helix found but nullptr returned!!\n"; continue; } HelixSeed tmp_helix_seed; @@ -305,15 +305,15 @@ namespace mu2e { if (helix_seed_vec.size() == 0) continue; -//----------------------------------------------------------------------------- -// now select the best helix to avoid duplicates -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // now select the best helix to avoid duplicates + //----------------------------------------------------------------------------- int index_best(-1); pickBestHelix(helix_seed_vec, index_best); -//----------------------------------------------------------------------------- -// fill seed information -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // fill seed information + //----------------------------------------------------------------------------- if ( (index_best>=0) && (index_best < 2) ){ Helicity hel_best = helix_seed_vec[index_best]._helix._helicity; if (_doSingleOutput) { @@ -338,9 +338,9 @@ namespace mu2e { // helix_seed_vec[index_best]._status.merge(TrkFitFlag::helixOK); // outseeds->push_back(helix_seed_vec[index_best]); if (_diagLevel > 0) { -//-------------------------------------------------------------------------------- -// fill diagnostic information -//-------------------------------------------------------------------------------- + //-------------------------------------------------------------------------------- + // fill diagnostic information + //-------------------------------------------------------------------------------- int nhitsMin(15); double mm2MeV = (3/10.)*_bz0; @@ -374,9 +374,9 @@ namespace mu2e { _data.shmeanr [loc] = _hfResult._diag.straw_mean_radius; _data.chi2d_helix [loc] = _hfResult._diag.chi2d_helix; if (_hfResult._diag.chi2d_helix>3) printf("[%s] : chi2Helix = %10.3f event number %8i\n", oname,_hfResult._diag.chi2d_helix,_iev); -//----------------------------------------------------------------------------- -// info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition -//----------------------------------------------------------------------------- + //----------------------------------------------------------------------------- + // info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition + //----------------------------------------------------------------------------- _data.loopId [loc] = _hfResult._diag.loopId_4; if (_hfResult._diag.loopId_4 == 1) { _data.chi2d_loop0 [loc] = _hfResult._diag.chi2_dof_circle_12; @@ -390,9 +390,9 @@ namespace mu2e { _data.npoints_loop1 [loc] = _hfResult._diag.n_active_11; } -//-------------------------------------------------------------------------------- -// info of the track candidate during the CAlHelixFinderAlg::findTrack loop -//-------------------------------------------------------------------------------- + //-------------------------------------------------------------------------------- + // info of the track candidate during the CAlHelixFinderAlg::findTrack loop + //-------------------------------------------------------------------------------- int counter(0); for (unsigned i=0; i<_hfResult._hitsUsed.size(); ++i){ if (_hfResult._hitsUsed[i] != 1) continue; @@ -419,22 +419,22 @@ namespace mu2e { // }//end panels loop // }//end faces loop } - else { - printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); + else { + printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); + } } - } - } -//-------------------------------------------------------------------------------- -// fill histograms -//-------------------------------------------------------------------------------- - if (_diagLevel > 0) { - _hmanager->fillHistograms(&_data); - } -//----------------------------------------------------------------------------- -// put reconstructed tracks into the event record -//----------------------------------------------------------------------------- - END:; + } + //-------------------------------------------------------------------------------- + // fill histograms + //-------------------------------------------------------------------------------- + if (_diagLevel > 0) { + _hmanager->fillHistograms(&_data); + } + //----------------------------------------------------------------------------- + // put reconstructed tracks into the event record + //----------------------------------------------------------------------------- +END:; int nseeds(0); if (_doSingleOutput) { nseeds += helcols[0]->size(); @@ -445,258 +445,258 @@ namespace mu2e { event.put(std::move(helcols[hel]),Helicity::name(hel)); } } - } - -//----------------------------------------------------------------------------- -// -//----------------------------------------------------------------------------- - void CalHelixFinder::endJob() { - // does this cause the file to close? - art::ServiceHandle tfs; - } - -//-------------------------------------------------------------------------------- -// set helix parameters -//----------------------------------------------------------------------------- - void CalHelixFinder::initHelixSeed(HelixSeed& HelSeed, CalHelixFinderData& HfResult) { - - HelixTraj* hel = HfResult.helix(); + } - double helixRadius = 1./fabs(hel->omega()); - double impactParam = hel->d0(); - double phi0 = hel->phi0(); - // double x0 = -(helixRadius + impactParam)*sin(phi0)*sig; - // double y0 = (helixRadius + impactParam)*cos(phi0)*sig; + //----------------------------------------------------------------------------- + // + //----------------------------------------------------------------------------- + void CalHelixFinder::endJob() { + // does this cause the file to close? + art::ServiceHandle tfs; + } - double x0 = -(1/hel->omega()+impactParam)*sin(phi0); - double y0 = (1/hel->omega()+impactParam)*cos(phi0); + //-------------------------------------------------------------------------------- + // set helix parameters + //----------------------------------------------------------------------------- + void CalHelixFinder::initHelixSeed(HelixSeed& HelSeed, CalHelixFinderData& HfResult) { - double dfdz = 1./hel->tanDip()/helixRadius; - double z0 = hel->z0(); - // center of the helix in the transverse plane - Hep3Vector center(x0, y0, 0); - //define the reconstructed helix parameters - HelSeed._helix._rcent = center.perp(); - HelSeed._helix._fcent = center.phi(); - HelSeed._helix._radius = helixRadius; - HelSeed._helix._lambda = 1./dfdz*_hfinder._dfdzsign; + HelixTraj* hel = HfResult.helix(); - HelSeed._helix._fz0 = phi0 - M_PI/2.*_hfinder._dfdzsign -z0*hel->omega()/hel->tanDip() ; + double helixRadius = 1./fabs(hel->omega()); + double impactParam = hel->d0(); + double phi0 = hel->phi0(); + // double x0 = -(helixRadius + impactParam)*sin(phi0)*sig; + // double y0 = (helixRadius + impactParam)*cos(phi0)*sig; - HelSeed._helix._helicity = HfResult._helicity; - HelSeed._status.merge(TrkFitFlag::CPRHelix); + double x0 = -(1/hel->omega()+impactParam)*sin(phi0); + double y0 = (1/hel->omega()+impactParam)*cos(phi0); - //include also the values of the chi2d - HelSeed._helix._chi2dXY = HfResult._sxy.chi2DofCircle(); - HelSeed._helix._chi2dZPhi = HfResult._szphi.chi2DofLine(); + double dfdz = 1./hel->tanDip()/helixRadius; + double z0 = hel->z0(); + // center of the helix in the transverse plane + Hep3Vector center(x0, y0, 0); + //define the reconstructed helix parameters + HelSeed._helix._rcent = center.perp(); + HelSeed._helix._fcent = center.phi(); + HelSeed._helix._radius = helixRadius; + HelSeed._helix._lambda = 1./dfdz*_hfinder._dfdzsign; - //now evaluate the helix T0 using the calorimeter cluster - double mm2MeV = (3/10.)*_bz0; - double tandip = hel->tanDip(); - double mom = helixRadius*mm2MeV/std::cos( std::atan(tandip)); - double beta = _tpart.beta(mom); + HelSeed._helix._fz0 = phi0 - M_PI/2.*_hfinder._dfdzsign -z0*hel->omega()/hel->tanDip() ; - double hel_t0; - double pitchAngle = M_PI/2. - std::atan(tandip); + HelSeed._helix._helicity = HfResult._helicity; + HelSeed._status.merge(TrkFitFlag::CPRHelix); - HelSeed._timeCluster = HfResult._timeClusterPtr; + //include also the values of the chi2d + HelSeed._helix._chi2dXY = HfResult._sxy.chi2DofCircle(); + HelSeed._helix._chi2dZPhi = HfResult._szphi.chi2DofLine(); - // cluster hits assigned to the reconsturcted Helix + //now evaluate the helix T0 using the calorimeter cluster + double mm2MeV = (3/10.)*_bz0; + double tandip = hel->tanDip(); + double mom = helixRadius*mm2MeV/std::cos( std::atan(tandip)); + double beta = _tpart.beta(mom); - int nhits = HfResult.nGoodHits(); - // printf("[CalHelixFinder::initHelixSeed] radius = %2.3f x0 = %2.3f y0 = %2.3f dfdz = %2.3e nhits = %i chi2XY = %2.3f chi2PHIZ = %2.3f\n", - // helixRadius, center.x(), center.y(), dfdz, nhits, HfResult._sxyw.chi2DofCircle(), HfResult._srphi.chi2DofLine()); - // printf("[CalHelixFinder::initHelixSeed] Index X Y Z PHI\n"); + double hel_t0; + double pitchAngle = M_PI/2. - std::atan(tandip); - // double z_start(0); - HelSeed._hhits.setParent(_chcol->parent()); + HelSeed._timeCluster = HfResult._timeClusterPtr; - for (int i=0; i_chHitsToProcess.at(hitInfo->panelHitIndex); - ComboHit hhit(*hit); - HelSeed._hhits.push_back(hhit); - } + // cluster hits assigned to the reconsturcted Helix - if(HfResult._timeClusterPtr->hasCaloCluster()) { - CLHEP::Hep3Vector gpos,tpos; + int nhits = HfResult.nGoodHits(); + // printf("[CalHelixFinder::initHelixSeed] radius = %2.3f x0 = %2.3f y0 = %2.3f dfdz = %2.3e nhits = %i chi2XY = %2.3f chi2PHIZ = %2.3f\n", + // helixRadius, center.x(), center.y(), dfdz, nhits, HfResult._sxyw.chi2DofCircle(), HfResult._srphi.chi2DofLine()); + // printf("[CalHelixFinder::initHelixSeed] Index X Y Z PHI\n"); - gpos = _hfinder._calorimeter->geomUtil().diskToMu2e(HfResult._timeClusterPtr->caloCluster()->diskID(), - HfResult._timeClusterPtr->caloCluster()->cog3Vector()); - tpos = _hfinder._calorimeter->geomUtil().mu2eToTracker(gpos); + // double z_start(0); + HelSeed._hhits.setParent(_chcol->parent()); - hel_t0 = HfResult._timeClusterPtr->caloCluster()->time() - (tpos.z() - z0)/sin(pitchAngle)/(beta*CLHEP::c_light); + for (int i=0; i_chHitsToProcess.at(hitInfo->panelHitIndex); + ComboHit hhit(*hit); + HelSeed._hhits.push_back(hhit); + } - } - else { + if(HfResult._timeClusterPtr->hasCaloCluster()) { + CLHEP::Hep3Vector gpos,tpos; - double tSum = 0.; + gpos = _hfinder._calorimeter->geomUtil().diskToMu2e(HfResult._timeClusterPtr->caloCluster()->diskID(), + HfResult._timeClusterPtr->caloCluster()->cog3Vector()); + tpos = _hfinder._calorimeter->geomUtil().mu2eToTracker(gpos); - size_t nHits = HelSeed.hits().size(); + hel_t0 = HfResult._timeClusterPtr->caloCluster()->time() - (tpos.z() - z0)/sin(pitchAngle)/(beta*CLHEP::c_light); - if(nHits == 0) { - hel_t0 = 0.; } else { - for(size_t i = 0; inhits(); - int ngoodhits(0); - // double minT(500.), maxT(2000.); - for (int i=0; ihits().at(i); - // StrawHitFlag flag = _shfcol->at(index); - ComboHit sh = _chcol ->at(index); - StrawHitFlag flag = sh.flag(); - int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); - if (bkg_hit) continue; - // if ( (sh.time() < minT) || (sh.time() > maxT) ) continue; - - ngoodhits += sh.nStrawHits(); - } + HelSeed._eDepAvg = HelSeed._hhits.eDepAvg(); - return ngoodhits; - } -//-------------------------------------------------------------------------------- -// function to select the best Helix among the results of the two helicity hypo -//-------------------------------------------------------------------------------- - void CalHelixFinder::pickBestHelix(std::vector& HelVec, int &Index_best){ - if (HelVec.size() == 1) { - Index_best = 0; - return; + //now set the HelixRecoDir + HelixTool ht(&HelSeed, _tracker); + float slope(0), slopeErr(0), chi2ndof(0); + ht.dirOfProp(slope, slopeErr, chi2ndof); + HelSeed._recoDir._slope = slope; + HelSeed._recoDir._slopeErr = slopeErr; + HelSeed._recoDir._chi2ndof = chi2ndof; } - const HelixSeed *h1, *h2; - const ComboHitCollection *tlist, *clist; - int nh1, nh2, natc(0); - const mu2e::HelixHit *hitt, *hitc; - - h1 = &HelVec[0]; -//------------------------------------------------------------------------------ -// check if an AlgorithmID collection has been created by the process -//----------------------------------------------------------------------------- - tlist = &h1->hits(); - nh1 = tlist->size(); - - h2 = &HelVec[1]; -//----------------------------------------------------------------------------- -// at Mu2e, 2 helices with different helicity could be duplicates of each other -//----------------------------------------------------------------------------- - clist = &h2->hits(); - nh2 = clist->size(); - -//----------------------------------------------------------------------------- -// check the number of common hits -//----------------------------------------------------------------------------- - for (int k=0; kat(k); - for (int l=0; lat(l); - if (hitt->index() == hitc->index()) { - natc += 1; - break; + //----------------------------------------------------------------------------- + int CalHelixFinder::initHelixFinderData(CalHelixFinderData& Data, + const TrkParticle& TPart, + const TrkFitDirection& FDir, + // const ComboHitCollection* ComboCollection , + // const StrawHitPositionCollection* ShPosCollection , + // const StrawHitFlagCollection* ShFlagCollection) { + const ComboHitCollection* ComboCollection ) { + Data._fit = TrkErrCode::fail; + Data._tpart = TPart; + Data._fdir = FDir; + + Data._chcol = ComboCollection; + // Data._shpos = ShPosCollection; + // Data._shfcol = ShFlagCollection; + + Data._radius = -1.0; + Data._dfdz = 0.; + Data._fz0 = 0.; + + return 0; } - } - } + int CalHelixFinder::goodHitsTimeCluster(const TimeCluster* TCluster){ + int nhits = TCluster->nhits(); + int ngoodhits(0); + // double minT(500.), maxT(2000.); + for (int i=0; ihits().at(i); + // StrawHitFlag flag = _shfcol->at(index); + ComboHit sh = _chcol ->at(index); + StrawHitFlag flag = sh.flag(); + int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); + if (bkg_hit) continue; + // if ( (sh.time() < minT) || (sh.time() > maxT) ) continue; + + ngoodhits += sh.nStrawHits(); + } - if ((natc > nh1/2.) || (natc > nh2/2.)) { + return ngoodhits; + } - //----------------------------------------------------------------------------- - // pick the helix with the largest number of hits - //----------------------------------------------------------------------------- - if (nh2 > nh1) { -//----------------------------------------------------------------------------- -// h2 is a winner, no need to save h1 -//----------------------------------------------------------------------------- - Index_best = 1; - return; - } - else if (nh1 > nh2){ -//----------------------------------------------------------------------------- -// h1 is a winner, mark h2 in hope that it will be OK, continue looping -//----------------------------------------------------------------------------- + //-------------------------------------------------------------------------------- + // function to select the best Helix among the results of the two helicity hypo + //-------------------------------------------------------------------------------- + void CalHelixFinder::pickBestHelix(std::vector& HelVec, int &Index_best){ + if (HelVec.size() == 1) { Index_best = 0; return; } -//----------------------------------------------------------------------------- -// in case they have the exact amount of hits, pick the one with better chi2dZphi -//----------------------------------------------------------------------------- - if (nh1 == nh2) { - float chi2dZphi_h1 = h1->helix().chi2dZPhi(); - float chi2dZphi_h2 = h2->helix().chi2dZPhi(); - if (chi2dZphi_h1 < chi2dZphi_h2){ - Index_best = 0; - return; - }else { + + const HelixSeed *h1, *h2; + const ComboHitCollection *tlist, *clist; + int nh1, nh2, natc(0); + const mu2e::HelixHit *hitt, *hitc; + + h1 = &HelVec[0]; + //------------------------------------------------------------------------------ + // check if an AlgorithmID collection has been created by the process + //----------------------------------------------------------------------------- + tlist = &h1->hits(); + nh1 = tlist->size(); + + h2 = &HelVec[1]; + //----------------------------------------------------------------------------- + // at Mu2e, 2 helices with different helicity could be duplicates of each other + //----------------------------------------------------------------------------- + clist = &h2->hits(); + nh2 = clist->size(); + + //----------------------------------------------------------------------------- + // check the number of common hits + //----------------------------------------------------------------------------- + for (int k=0; kat(k); + for (int l=0; lat(l); + if (hitt->index() == hitc->index()) { + natc += 1; + break; + } + } + } + + + if ((natc > nh1/2.) || (natc > nh2/2.)) { + + //----------------------------------------------------------------------------- + // pick the helix with the largest number of hits + //----------------------------------------------------------------------------- + if (nh2 > nh1) { + //----------------------------------------------------------------------------- + // h2 is a winner, no need to save h1 + //----------------------------------------------------------------------------- Index_best = 1; return; } + else if (nh1 > nh2){ + //----------------------------------------------------------------------------- + // h1 is a winner, mark h2 in hope that it will be OK, continue looping + //----------------------------------------------------------------------------- + Index_best = 0; + return; + } + //----------------------------------------------------------------------------- + // in case they have the exact amount of hits, pick the one with better chi2dZphi + //----------------------------------------------------------------------------- + if (nh1 == nh2) { + float chi2dZphi_h1 = h1->helix().chi2dZPhi(); + float chi2dZphi_h2 = h2->helix().chi2dZPhi(); + if (chi2dZphi_h1 < chi2dZphi_h2){ + Index_best = 0; + return; + }else { + Index_best = 1; + return; + } + } + }else { + //----------------------------------------------------------------------------- + // this is the case where we consider the two helices independent, so we want + // to store both + //----------------------------------------------------------------------------- + Index_best = 2; + return; } - }else { -//----------------------------------------------------------------------------- -// this is the case where we consider the two helices independent, so we want -// to store both -//----------------------------------------------------------------------------- - Index_best = 2; - return; - } - } + } -} + } -using mu2e::CalHelixFinder; -DEFINE_ART_MODULE(CalHelixFinder) + using mu2e::CalHelixFinder; + DEFINE_ART_MODULE(CalHelixFinder) diff --git a/GeometryService/inc/GeometryService.hh b/GeometryService/inc/GeometryService.hh index 806c80d2e2..f184ec1b02 100644 --- a/GeometryService/inc/GeometryService.hh +++ b/GeometryService/inc/GeometryService.hh @@ -21,6 +21,7 @@ #include "fhiclcpp/types/Table.h" #include "cetlib_except/exception.h" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/ConfigTools/inc/SimpleConfig.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" #include "boost/shared_ptr.hpp" @@ -48,6 +49,7 @@ public: fhicl::Atom tool_type{Name("tool_type"),"Mu2e"}; }; + using KKMaterialConfig = KKMaterial::Config; struct Config { using Name=fhicl::Name; using Comment=fhicl::Comment; @@ -60,6 +62,7 @@ public: fhicl::Atom printConfig{Name("printConfig"),false}; fhicl::Atom printConfigTopLevel{Name("printConfigTopLevel"),false}; fhicl::Table simulatedDetector{Name("simulatedDetector")}; + fhicl::Table matSettings{Name("KinKalMaterial")}; }; using Parameters= art::ServiceTable; @@ -129,6 +132,7 @@ private: std::unique_ptr _bfConfig; const fhicl::ParameterSet _simulatedDetector; + const KKMaterialConfig _kkMat; // Load G4 geometry options std::unique_ptr _g4GeomOptions; diff --git a/GeometryService/inc/KinKalGeomMaker.hh b/GeometryService/inc/KinKalGeomMaker.hh index d0adbe8dbc..9c89168c9b 100644 --- a/GeometryService/inc/KinKalGeomMaker.hh +++ b/GeometryService/inc/KinKalGeomMaker.hh @@ -5,6 +5,7 @@ // Original author: Dave Brown (LBNL) 4/2026 // #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" namespace mu2e { class KinKalGeomMaker { public: diff --git a/GeometryService/src/GeometryService.cc b/GeometryService/src/GeometryService.cc index 67500f6694..f491783b12 100644 --- a/GeometryService/src/GeometryService.cc +++ b/GeometryService/src/GeometryService.cc @@ -61,6 +61,7 @@ #include "Offline/BeamlineGeom/inc/TSdA.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/GeometryService/inc/DiskCalorimeterMaker.hh" #include "Offline/CalorimeterGeom/inc/DiskCalorimeter.hh" @@ -96,7 +97,7 @@ using namespace std; namespace mu2e { GeometryService::GeometryService( const Parameters& pars, - art::ActivityRegistry&iRegistry) : + art::ActivityRegistry&iRegistry) : _inputfile( pars().inputFile()), _bFieldFile( pars().bFieldFile()), _allowReplacement( pars().allowReplacement()), @@ -107,6 +108,7 @@ namespace mu2e { _printTopLevel( pars().printConfigTopLevel()), _config(nullptr), _simulatedDetector( pars.get_PSet().get("simulatedDetector")), + _kkMat( pars().matSettings()), _standardMu2eDetector( _simulatedDetector.get("tool_type") == "Mu2e"), _detectors() { @@ -120,246 +122,244 @@ namespace mu2e { // This template can be defined here because this is a private method which is only // used by the code below in the same file. template - void GeometryService::addDetector(std::unique_ptr d) - { - if(_detectors.find(typeid(DET).name())!=_detectors.end()) { - throw cet::exception("GEOM") << "failed to install detector with type name " - << typeid(DET).name() << "\n"; - } + void GeometryService::addDetector(std::unique_ptr d) + { + if(_detectors.find(typeid(DET).name())!=_detectors.end()) { + throw cet::exception("GEOM") << "failed to install detector with type name " + << typeid(DET).name() << "\n"; + } DetectorPtr ptr(d.release()); _detectors[typeid(DET).name()] = ptr; - } + } template - void GeometryService::addDetectorAliasToBaseClass(std::unique_ptr d) - { + void GeometryService::addDetectorAliasToBaseClass(std::unique_ptr d) + { - std::string OriginalName = typeid(DET).name(); - DetMap::iterator it(_detectors.find(OriginalName)); + std::string OriginalName = typeid(DET).name(); + DetMap::iterator it(_detectors.find(OriginalName)); - if(it==_detectors.end()) - throw cet::exception("GEOM") - << "Can not alias an inexistant detector, detector " << OriginalName << "\n"; + if(it==_detectors.end()) + throw cet::exception("GEOM") + << "Can not alias an inexistant detector, detector " << OriginalName << "\n"; - std::string detectorName= typeid(DETALIAS).name() ; - _detectors[detectorName] = it->second; - } + std::string detectorName= typeid(DETALIAS).name() ; + _detectors[detectorName] = it->second; + } void - GeometryService::preBeginRun(art::Run const &) { - - if(++_run_count > 1) { - return; - } + GeometryService::preBeginRun(art::Run const &) { - _config = unique_ptr(new SimpleConfig(_inputfile, - _allowReplacement, - _messageOnReplacement, - _messageOnDefault )); - _config->printOpen(cout,"Geometry"); - - _bfConfig = unique_ptr(new SimpleConfig(_bFieldFile, - _allowReplacement, - _messageOnReplacement, - _messageOnDefault )); - _bfConfig->printOpen(cout,"BField"); - - - if(_printTopLevel) { - //print the top level geometry file contents - //the top level often contains a single named config file or a list of specific version files - ConfigFileLookupPolicy configFile; - std::string file = configFile(_inputfile); - std::ifstream in(file.c_str()); - if ( !in ) { - // No conf file for this test. - throw cet::exception("Geom") - << "GeometryService: Cannot open input file: " - << file - << endl; + if(++_run_count > 1) { + return; } - std::cout << "GeometryService: printing top level geometry file:\n"; - std::string line; - while ( in ){ - std::getline(in,line); - if ( !in ){ - break; - } - std::cout << line.c_str() << std::endl; + _config = unique_ptr(new SimpleConfig(_inputfile, + _allowReplacement, + _messageOnReplacement, + _messageOnDefault )); + _config->printOpen(cout,"Geometry"); + + _bfConfig = unique_ptr(new SimpleConfig(_bFieldFile, + _allowReplacement, + _messageOnReplacement, + _messageOnDefault )); + _bfConfig->printOpen(cout,"BField"); + + + if(_printTopLevel) { + //print the top level geometry file contents + //the top level often contains a single named config file or a list of specific version files + ConfigFileLookupPolicy configFile; + std::string file = configFile(_inputfile); + std::ifstream in(file.c_str()); + if ( !in ) { + // No conf file for this test. + throw cet::exception("Geom") + << "GeometryService: Cannot open input file: " + << file + << endl; + } + std::cout << "GeometryService: printing top level geometry file:\n"; + std::string line; + while ( in ){ + std::getline(in,line); + if ( !in ){ + break; + } + + std::cout << line.c_str() << std::endl; + } + std::cout << "GeometryService: finished printing top level geometry file.\n"; } - std::cout << "GeometryService: finished printing top level geometry file.\n"; - } - // Print final state of file after all substitutions. - if ( _printConfig ){ _config->print(cout, "Geom: "); } - - // 2019-03-24 P.M. : *not needed* decide if this is standard Mu2e detector or something else ... - - if (!isStandardMu2eDetector() || - !_config->getBool("mu2e.standardDetector",true)) { - cout << "Non-standard mu2e configuration, assuming it is intentional" << endl; - return; - } - - // Initialize geometry options - _g4GeomOptions = unique_ptr( new G4GeometryOptions( *_config ) ); + // Print final state of file after all substitutions. + if ( _printConfig ){ _config->print(cout, "Geom: "); } - // Throw if the configuration is not self consistent. - checkConfig(); + // 2019-03-24 P.M. : *not needed* decide if this is standard Mu2e detector or something else ... - // This must be the first detector added since other makers may wish to use it. - std::unique_ptr tmpDetSys(DetectorSystemMaker::make(*_config)); - const DetectorSystem& detSys = *tmpDetSys.get(); - addDetector(std::move(tmpDetSys)); - - // Make a detector for every component present in the configuration. - - std::unique_ptr tmpBeamline(BeamlineMaker::make(*_config)); - const Beamline& beamline = *tmpBeamline.get(); - addDetector(std::move(tmpBeamline)); + if (!isStandardMu2eDetector() || + !_config->getBool("mu2e.standardDetector",true)) { + cout << "Non-standard mu2e configuration, assuming it is intentional" << endl; + return; + } - std::unique_ptr tmpProdTgt(ProductionTargetMaker::make(*_config, beamline.solenoidOffset())); - const ProductionTarget& prodTarget = *tmpProdTgt.get(); - addDetector(std::move(tmpProdTgt)); + // Initialize geometry options + _g4GeomOptions = unique_ptr( new G4GeometryOptions( *_config ) ); - std::unique_ptr - tmpProductionSolenoid(ProductionSolenoidMaker(*_config, beamline.solenoidOffset()).getProductionSolenoidPtr()); + // Throw if the configuration is not self consistent. + checkConfig(); - const ProductionSolenoid& ps = *tmpProductionSolenoid.get(); - addDetector(std::move(tmpProductionSolenoid)); + // This must be the first detector added since other makers may wish to use it. + std::unique_ptr tmpDetSys(DetectorSystemMaker::make(*_config)); + const DetectorSystem& detSys = *tmpDetSys.get(); + addDetector(std::move(tmpDetSys)); - std::unique_ptr - tmpPSE(PSEnclosureMaker::make(*_config, ps.psEndRefPoint())); - const PSEnclosure& pse = *tmpPSE.get(); - addDetector(std::move(tmpPSE)); + // Make a detector for every component present in the configuration. - // The Z coordinate of the boundary between PS and TS vacua - StraightSection const * ts1vac = beamline.getTS().getTSVacuum( TransportSolenoid::TSRegion::TS1 ); - const double vacPS_TS_z = ts1vac->getGlobal().z() - ts1vac->getHalfLength(); + std::unique_ptr tmpBeamline(BeamlineMaker::make(*_config)); + const Beamline& beamline = *tmpBeamline.get(); + addDetector(std::move(tmpBeamline)); - addDetector(PSVacuumMaker::make(*_config, ps, pse, vacPS_TS_z)); + std::unique_ptr tmpProdTgt(ProductionTargetMaker::make(*_config, beamline.solenoidOffset())); + const ProductionTarget& prodTarget = *tmpProdTgt.get(); + addDetector(std::move(tmpProdTgt)); - // Use ProductionTarget's built-in method to get the correct position based on model type - addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.targetPositionByVersion())); + std::unique_ptr + tmpProductionSolenoid(ProductionSolenoidMaker(*_config, beamline.solenoidOffset()).getProductionSolenoidPtr()); + const ProductionSolenoid& ps = *tmpProductionSolenoid.get(); + addDetector(std::move(tmpProductionSolenoid)); + std::unique_ptr + tmpPSE(PSEnclosureMaker::make(*_config, ps.psEndRefPoint())); + const PSEnclosure& pse = *tmpPSE.get(); + addDetector(std::move(tmpPSE)); + // The Z coordinate of the boundary between PS and TS vacua + StraightSection const * ts1vac = beamline.getTS().getTSVacuum( TransportSolenoid::TSRegion::TS1 ); + const double vacPS_TS_z = ts1vac->getGlobal().z() - ts1vac->getHalfLength(); + addDetector(PSVacuumMaker::make(*_config, ps, pse, vacPS_TS_z)); + // Use ProductionTarget's built-in method to get the correct position based on model type + addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.targetPositionByVersion())); - // Construct building solids - std::unique_ptr tmphall(Mu2eHallMaker::makeBuilding(*_g4GeomOptions,*_config)); - const Mu2eHall& hall = *tmphall.get(); + // Construct building solids + std::unique_ptr tmphall(Mu2eHallMaker::makeBuilding(*_g4GeomOptions,*_config)); + const Mu2eHall& hall = *tmphall.get(); - // Determine Mu2e envelope from building solids - std::unique_ptr mu2eEnv (new Mu2eEnvelope(hall,*_config)); + // Determine Mu2e envelope from building solids + std::unique_ptr mu2eEnv (new Mu2eEnvelope(hall,*_config)); - // Make dirt based on Mu2e envelope - Mu2eHallMaker::makeDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - Mu2eHallMaker::makeRotated( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - Mu2eHallMaker::makeTrapDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + // Make dirt based on Mu2e envelope + Mu2eHallMaker::makeDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + Mu2eHallMaker::makeRotated( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + Mu2eHallMaker::makeTrapDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - addDetector(std::move( tmphall ) ); - addDetector(std::move( mu2eEnv ) ); + addDetector(std::move( tmphall ) ); + addDetector(std::move( mu2eEnv ) ); - std::unique_ptr tmpDump(ProtonBeamDumpMaker::make(*_config, hall)); - const ProtonBeamDump& dump = *tmpDump.get(); - addDetector(std::move(tmpDump)); + std::unique_ptr tmpDump(ProtonBeamDumpMaker::make(*_config, hall)); + const ProtonBeamDump& dump = *tmpDump.get(); + addDetector(std::move(tmpDump)); - // beamline info used to position DS - std::unique_ptr tmpDS( DetectorSolenoidMaker::make( *_config, beamline ) ); - const DetectorSolenoid& ds = *tmpDS.get(); - addDetector(std::move(tmpDS)); + // beamline info used to position DS + std::unique_ptr tmpDS( DetectorSolenoidMaker::make( *_config, beamline ) ); + const DetectorSolenoid& ds = *tmpDS.get(); + addDetector(std::move(tmpDS)); - // DS info used to position DS downstream shielding - addDetector( DetectorSolenoidShieldingMaker::make( *_config, ds ) ); + // DS info used to position DS downstream shielding + addDetector( DetectorSolenoidShieldingMaker::make( *_config, ds ) ); - std::unique_ptr tmptgt(StoppingTargetMaker(detSys.getOrigin(), *_config).getTargetPtr()); - const StoppingTarget& target = *tmptgt.get(); - addDetector(std::move(tmptgt)); + std::unique_ptr tmptgt(StoppingTargetMaker(detSys.getOrigin(), *_config).getTargetPtr()); + const StoppingTarget& target = *tmptgt.get(); + addDetector(std::move(tmptgt)); - if (_config->getBool("hasTracker",false)){ - TrackerMaker ttm( *_config ); - addDetector( ttm.getTrackerPtr() ); - } + if (_config->getBool("hasTracker",false)){ + TrackerMaker ttm( *_config ); + addDetector( ttm.getTrackerPtr() ); + } - if(_config->getBool("hasMBS",false)){ - MBSMaker mbs( *_config, beamline.solenoidOffset() ); - addDetector( mbs.getMBSPtr() ); - } + if(_config->getBool("hasMBS",false)){ + MBSMaker mbs( *_config, beamline.solenoidOffset() ); + addDetector( mbs.getMBSPtr() ); + } - if(_config->getBool("hasDiskCalorimeter",false)){ - DiskCalorimeterMaker calorm( *_config, beamline.solenoidOffset() ); - addDetector( calorm.calorimeterPtr() ); - addDetectorAliasToBaseClass( calorm.calorimeterPtr() ); //add an alias to detector list - } + if(_config->getBool("hasDiskCalorimeter",false)){ + DiskCalorimeterMaker calorm( *_config, beamline.solenoidOffset() ); + addDetector( calorm.calorimeterPtr() ); + addDetectorAliasToBaseClass( calorm.calorimeterPtr() ); //add an alias to detector list + } - if(_config->getBool("hasCosmicRayShield",false)){ - CosmicRayShieldMaker crs( *_config, beamline.solenoidOffset() ); - addDetector( crs.getCosmicRayShieldPtr() ); - } + if(_config->getBool("hasCosmicRayShield",false)){ + CosmicRayShieldMaker crs( *_config, beamline.solenoidOffset() ); + addDetector( crs.getCosmicRayShieldPtr() ); + } - if(_config->getBool("hasTSdA",false)){ - addDetector( TSdAMaker::make(*_config,ds) ); - } + if(_config->getBool("hasTSdA",false)){ + addDetector( TSdAMaker::make(*_config,ds) ); + } - if(_config->getBool("hasExternalShielding",false)) { - addDetector( ExtShieldUpstreamMaker::make(*_config) ); - addDetector( ExtShieldDownstreamMaker::make(*_config)); - addDetector( SaddleMaker::make(*_config)); - addDetector( PipeMaker::make(*_config)); - addDetector( ElectronicRackMaker::make(*_config)); - } + if(_config->getBool("hasExternalShielding",false)) { + addDetector( ExtShieldUpstreamMaker::make(*_config) ); + addDetector( ExtShieldDownstreamMaker::make(*_config)); + addDetector( SaddleMaker::make(*_config)); + addDetector( PipeMaker::make(*_config)); + addDetector( ElectronicRackMaker::make(*_config)); + } - std::unique_ptr tmpemb(ExtMonFNALBuildingMaker::make(*_config, hall, dump)); - const ExtMonFNALBuilding& emfb = *tmpemb.get(); - addDetector(std::move(tmpemb)); - if(_config->getBool("hasExtMonFNAL",false)){ - addDetector(ExtMonFNAL::ExtMonMaker::make(*_config, emfb)); - } + std::unique_ptr tmpemb(ExtMonFNALBuildingMaker::make(*_config, hall, dump)); + const ExtMonFNALBuilding& emfb = *tmpemb.get(); + addDetector(std::move(tmpemb)); + if(_config->getBool("hasExtMonFNAL",false)){ + addDetector(ExtMonFNAL::ExtMonMaker::make(*_config, emfb)); + } - if (_config->getBool("hasPTM",false) ){ - std::unique_ptr ptmon(PTMMaker::make(*_config)); - addDetector(std::move(ptmon)); - } + if (_config->getBool("hasPTM",false) ){ + std::unique_ptr ptmon(PTMMaker::make(*_config)); + addDetector(std::move(ptmon)); + } - if(_config->getBool("hasSTM",false)){ - STMMaker stm( *_config, beamline.solenoidOffset() ); - addDetector( stm.getSTMPtr() ); - } + if(_config->getBool("hasSTM",false)){ + STMMaker stm( *_config, beamline.solenoidOffset() ); + addDetector( stm.getSTMPtr() ); + } - if(_config->getBool("hasVirtualDetector",false)){ - addDetector(VirtualDetectorMaker::make(*_config)); - } + if(_config->getBool("hasVirtualDetector",false)){ + addDetector(VirtualDetectorMaker::make(*_config)); + } - if(_bfConfig->getBool("hasBFieldManager",false)){ - std::unique_ptr bfc( BFieldConfigMaker(*_bfConfig, beamline).getBFieldConfig() ); - BFieldManagerMaker bfmgr(*bfc); - addDetector(std::move(bfc)); - addDetector(bfmgr.getBFieldManager()); - } + if(_bfConfig->getBool("hasBFieldManager",false)){ + std::unique_ptr bfc( BFieldConfigMaker(*_bfConfig, beamline).getBFieldConfig() ); + BFieldManagerMaker bfmgr(*bfc); + addDetector(std::move(bfc)); + addDetector(bfmgr.getBFieldManager()); + } - if(_config->getBool("hasProtonAbsorber",false) && !_config->getBool("protonabsorber.isHelical", false) ){ - MECOStyleProtonAbsorberMaker mecopam( *_config, ds, target); - addDetector( mecopam.getMECOStyleProtonAbsorberPtr() ); - } + if(_config->getBool("hasProtonAbsorber",false) && !_config->getBool("protonabsorber.isHelical", false) ){ + MECOStyleProtonAbsorberMaker mecopam( *_config, ds, target); + addDetector( mecopam.getMECOStyleProtonAbsorberPtr() ); + } - // This class has a default c'tor with all available information internally. - std::unique_ptr dusafMu2e{ std::make_unique() }; - addDetector( std::move(dusafMu2e) ); + // This class has a default c'tor with all available information internally. + std::unique_ptr dusafMu2e{ std::make_unique() }; + addDetector( std::move(dusafMu2e) ); - // build KinKalGeom, used in track reconstruction and extrapolation - KinKalGeomMaker kkgm; - addDetector( std::move(kkgm.makeKKG()) ); + // build KinKalGeom, used in track reconstruction and extrapolation + KinKalGeomMaker kkgm; + addDetector( std::move(kkgm.makeKKG()) ); + // directly build KKMaterial; it's constructor does everything + addDetector( std::make_unique(_kkMat)); + return ; - } // preBeginRun() + } // preBeginRun() // Check that the configuration is self consistent. void GeometryService::checkConfig(){ diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index 5c955ef72e..683b22a997 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -79,11 +79,41 @@ namespace mu2e { } void KinKalGeomMaker::makeDS() { - // currently use hard-coded geometry - auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),950,5450); - auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),1328,5450); // bounding surfaces + GeomHandle det; + GeomHandle ds; +// std::cout << "DS Cryo or " << ds->rOut1() << "," << ds->rOut2() << " ir " << ds->rIn1()<<","<< ds->rIn2() << " halfl " << ds->halfLength() +// << " zpos " << ds->position().z() << " material " << ds->material() << std::endl; +// std::cout << "DS shield or " << ds->shield_rOut1() << "," << ds->shield_rOut2() << " ir " << ds->shield_rIn1()<<","<< ds->shield_rIn2() << " halfl " << ds->shield_halfLength() << " material " << ds->shield_material() << std::endl; +// std::cout << "DS ncoils " << ds->nCoils() << std::endl; +// for(size_t icoil = 0; icoil < static_cast(ds->nCoils()); icoil++){ +// std::cout << "DS coil ir " << ds->coil_rIn() << " or " << ds->coil_rOut()[icoil] << " length " << ds->coil_zLength()[icoil] << " zpos " << ds->coil_zPosition()[icoil] +// << " material " << ds->coil_materials()[icoil] << std::endl; +// } + //DS Cryo or 1303,1328 ir 950,969.05 halfl 5450 zpos 8689 material StainlessSteel + //DS shield or 1237.3,1250 ir 1010,1022.7 halfl 5287.7 material G4_Al + //DS ncoils 11 + //DS coil ir 1050 or 1091 length 419.75 zpos 3539 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 3964 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 4389 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 5042 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 5699 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 6369 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 7176 material DS1CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 7949 material DS2CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 9761 material DS2CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 11544 material DS2CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 13326 material DS1CoilMix + // + // + // + // + auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),ds->rIn1(),ds->halfLength()); + auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),ds->rOut2(),ds->halfLength()); // bounding surfaces auto front= std::make_shared(outer->frontDisk()); auto back= std::make_shared(outer->backDisk()); + + + // hard-coded for now auto ipa= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-2770),300.0,500.0); auto ipafront= std::make_shared(ipa->frontDisk()); auto ipaback= std::make_shared(ipa->backDisk()); diff --git a/GeometryService/src/ProductionTargetMaker.cc b/GeometryService/src/ProductionTargetMaker.cc index fbc1161300..5d053843a9 100644 --- a/GeometryService/src/ProductionTargetMaker.cc +++ b/GeometryService/src/ProductionTargetMaker.cc @@ -21,13 +21,13 @@ namespace mu2e { std::unique_ptr ProductionTargetMaker::make(const SimpleConfig& c, double solenoidOffset) { const std::string modelName = c.getString("targetPS_model"); const auto& registry = getMakerRegistry(); - + auto it = registry.find(modelName); if (it != registry.end()) { return it->second(c, solenoidOffset); } - - throw cet::exception("GEOM") + + throw cet::exception("GEOM") << "illegal production target model specified = " << modelName << std::endl; } @@ -41,9 +41,9 @@ namespace mu2e { c.getDouble("targetPS_rotX") * CLHEP::degree, c.getDouble("targetPS_rotY") * CLHEP::degree, CLHEP::Hep3Vector(solenoidOffset, - 0, - c.getDouble("productionTarget.zNominal") - ) + 0, + c.getDouble("productionTarget.zNominal") + ) + c.getHep3Vector("productionTarget.offset", CLHEP::Hep3Vector(0,0,0)), c.getInt ("targetPS_nFins" , 0 ), @@ -55,7 +55,7 @@ namespace mu2e { c.getDouble("targetPS_hubAngleDS" ,0.0), c.getDouble("targetPS_hubOverhangUS",0.0), c.getDouble("targetPS_hubOverhangDS",0.0) ) ); - double trgtMaxAngle = c.getDouble("targetPS_rotY"); + double trgtMaxAngle = c.getDouble("targetPS_rotY"); if (c.getDouble("targetPS_rotX")>trgtMaxAngle) { trgtMaxAngle=c.getDouble("targetPS_rotX"); } trgtMaxAngle *= CLHEP::deg; @@ -79,196 +79,196 @@ namespace mu2e { //double Hub_hang_TotLength = Hub_hang_Length; if (Hub_overhang_angle>=0.0 && Hub_overhang_angle<90.0) { - Hub_overhang_angle *= CLHEP::deg; - double deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); - double cosHub_overhang_angle = cos(Hub_overhang_angle); - double appThick = Hub_thickness/cosHub_overhang_angle; - double actual_hang_Length = Hub_hang_Length; - if ( version > 1 ) actual_hang_Length = Hub_thickness/tan(Hub_overhang_angle); - tgtPS->setHubLenDS(actual_hang_Length); - //Hub_hang_TotLength+=Hub_overhang_Length; - HubRgtCornersZ.push_back(tgtPS->halfLength()+Hub_overhang_Length - - tgtPS->hubDistDS()); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()+deltaRad); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+appThick); - - HubRgtCornersZ.push_back(tgtPS->halfLength()-tgtPS->hubDistDS()); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); - - - HubRgtCornersZ.push_back(tgtPS->halfLength() - actual_hang_Length - - tgtPS->hubDistDS()); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - - if ( version > 1 ) { - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+0.001); - } else { - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); - } - // Finished downstream/right, now do upstream/left - if (version > 1) { - Hub_overhang_Length = tgtPS->hubOverhangUS(); - Hub_overhang_angle = tgtPS->hubAngleUS(); - Hub_overhang_angle *= CLHEP::deg; - deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); - cosHub_overhang_angle = cos(Hub_overhang_angle); - appThick = Hub_thickness/cosHub_overhang_angle; - actual_hang_Length = Hub_hang_Length; - if ( version > 1 ) actual_hang_Length = Hub_thickness/tan(Hub_overhang_angle); - tgtPS->setHubLenUS(actual_hang_Length); - } - - HubLftCornersZ.push_back(-tgtPS->halfLength()-Hub_overhang_Length - + tgtPS->hubDistUS()); - HubLftCornersInnRadii.push_back(tgtPS->rOut()+deltaRad); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+appThick); - - HubLftCornersZ.push_back(-tgtPS->halfLength() + tgtPS->hubDistUS()); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); - - HubLftCornersZ.push_back(-tgtPS->halfLength()+actual_hang_Length - + tgtPS->hubDistUS()); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - if (version > 1 ) { - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+0.1); - } else { - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); - } - - // Created both hubs - - - double SpkAnchrPosZ = HubRgtCornersZ.at(0)-spokeAnchordist*cosHub_overhang_angle; - double SpkAnchrPosRad = HubRgtCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); - double iSpkAngleRgt(0.0), iSpkAngleLft(-spokeSideDangle), tmpAngle(0.0); - for (int iSpk=0; iSpk 1) { - Hub_overhang_Length = tgtPS->hubOverhangDS(); - Hub_overhang_angle = tgtPS->hubAngleDS(); - Hub_overhang_angle *= CLHEP::deg; - deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); - cosHub_overhang_angle = cos(Hub_overhang_angle); - appThick = Hub_thickness/cosHub_overhang_angle; - SpkAnchrPosZ = HubRgtCornersZ.at(0)-spokeAnchordist*cosHub_overhang_angle; - SpkAnchrPosRad = HubRgtCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); - } - tmpAngle = iSpkAngleRgt * CLHEP::deg; - std::map::iterator pntPos; - pntPos = tgtPS->_anchoringPntsRgt.insert( - std::pair( - iSpkAngleRgt, - CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), - SpkAnchrPosRad*sin(tmpAngle), - SpkAnchrPosZ) - ) - ).first; - pntPos->second.transform(tgtPS->protonBeamRotation()); - // Switch back to Upstream version of #s - // if version 2 or higher... - if (version > 1) { - Hub_overhang_Length = tgtPS->hubOverhangUS(); - Hub_overhang_angle = tgtPS->hubAngleUS(); - Hub_overhang_angle *= CLHEP::deg; - deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); - cosHub_overhang_angle = cos(Hub_overhang_angle); - appThick = Hub_thickness/cosHub_overhang_angle; - SpkAnchrPosZ = HubLftCornersZ.at(0)+spokeAnchordist*cosHub_overhang_angle; - SpkAnchrPosZ *= -1.0; - SpkAnchrPosRad = HubLftCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); - } - tmpAngle = iSpkAngleLft * CLHEP::deg; - pntPos = tgtPS->_anchoringPntsLft.insert( - std::pair( - iSpkAngleLft, - CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), - SpkAnchrPosRad*sin(tmpAngle), - -SpkAnchrPosZ) - ) - ).first; - pntPos->second.transform(tgtPS->protonBeamRotation()); - iSpkAngleRgt += spokeAngleStep; - iSpkAngleLft += spokeAngleStep; - } - - if ( version > 1 ) { - double maxHang = tgtPS->hubOverhangUS(); - double tmpHang = tgtPS->hubOverhangDS(); - if ( tmpHang > maxHang ) maxHang = tmpHang; - tgtPS->_envelHalfLength += maxHang + 2.0; - } else { - tgtPS->_envelHalfLength += Hub_overhang_Length; - } + Hub_overhang_angle *= CLHEP::deg; + double deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); + double cosHub_overhang_angle = cos(Hub_overhang_angle); + double appThick = Hub_thickness/cosHub_overhang_angle; + double actual_hang_Length = Hub_hang_Length; + if ( version > 1 ) actual_hang_Length = Hub_thickness/tan(Hub_overhang_angle); + tgtPS->setHubLenDS(actual_hang_Length); + //Hub_hang_TotLength+=Hub_overhang_Length; + HubRgtCornersZ.push_back(tgtPS->halfLength()+Hub_overhang_Length + - tgtPS->hubDistDS()); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()+deltaRad); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+appThick); + + HubRgtCornersZ.push_back(tgtPS->halfLength()-tgtPS->hubDistDS()); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); + + + HubRgtCornersZ.push_back(tgtPS->halfLength() - actual_hang_Length + - tgtPS->hubDistDS()); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + + if ( version > 1 ) { + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+0.001); + } else { + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); + } + // Finished downstream/right, now do upstream/left + if (version > 1) { + Hub_overhang_Length = tgtPS->hubOverhangUS(); + Hub_overhang_angle = tgtPS->hubAngleUS(); + Hub_overhang_angle *= CLHEP::deg; + deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); + cosHub_overhang_angle = cos(Hub_overhang_angle); + appThick = Hub_thickness/cosHub_overhang_angle; + actual_hang_Length = Hub_hang_Length; + if ( version > 1 ) actual_hang_Length = Hub_thickness/tan(Hub_overhang_angle); + tgtPS->setHubLenUS(actual_hang_Length); + } + + HubLftCornersZ.push_back(-tgtPS->halfLength()-Hub_overhang_Length + + tgtPS->hubDistUS()); + HubLftCornersInnRadii.push_back(tgtPS->rOut()+deltaRad); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+appThick); + + HubLftCornersZ.push_back(-tgtPS->halfLength() + tgtPS->hubDistUS()); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); + + HubLftCornersZ.push_back(-tgtPS->halfLength()+actual_hang_Length + + tgtPS->hubDistUS()); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + if (version > 1 ) { + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+0.1); + } else { + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); + } + + // Created both hubs + + + double SpkAnchrPosZ = HubRgtCornersZ.at(0)-spokeAnchordist*cosHub_overhang_angle; + double SpkAnchrPosRad = HubRgtCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); + double iSpkAngleRgt(0.0), iSpkAngleLft(-spokeSideDangle), tmpAngle(0.0); + for (int iSpk=0; iSpk 1) { + Hub_overhang_Length = tgtPS->hubOverhangDS(); + Hub_overhang_angle = tgtPS->hubAngleDS(); + Hub_overhang_angle *= CLHEP::deg; + deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); + cosHub_overhang_angle = cos(Hub_overhang_angle); + appThick = Hub_thickness/cosHub_overhang_angle; + SpkAnchrPosZ = HubRgtCornersZ.at(0)-spokeAnchordist*cosHub_overhang_angle; + SpkAnchrPosRad = HubRgtCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); + } + tmpAngle = iSpkAngleRgt * CLHEP::deg; + std::map::iterator pntPos; + pntPos = tgtPS->_anchoringPntsRgt.insert( + std::pair( + iSpkAngleRgt, + CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), + SpkAnchrPosRad*sin(tmpAngle), + SpkAnchrPosZ) + ) + ).first; + pntPos->second.transform(tgtPS->protonBeamRotation()); + // Switch back to Upstream version of #s + // if version 2 or higher... + if (version > 1) { + Hub_overhang_Length = tgtPS->hubOverhangUS(); + Hub_overhang_angle = tgtPS->hubAngleUS(); + Hub_overhang_angle *= CLHEP::deg; + deltaRad = Hub_overhang_Length*tan(Hub_overhang_angle); + cosHub_overhang_angle = cos(Hub_overhang_angle); + appThick = Hub_thickness/cosHub_overhang_angle; + SpkAnchrPosZ = HubLftCornersZ.at(0)+spokeAnchordist*cosHub_overhang_angle; + SpkAnchrPosZ *= -1.0; + SpkAnchrPosRad = HubLftCornersOutRadii.at(0)-spokeAnchordist*sin(Hub_overhang_angle); + } + tmpAngle = iSpkAngleLft * CLHEP::deg; + pntPos = tgtPS->_anchoringPntsLft.insert( + std::pair( + iSpkAngleLft, + CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), + SpkAnchrPosRad*sin(tmpAngle), + -SpkAnchrPosZ) + ) + ).first; + pntPos->second.transform(tgtPS->protonBeamRotation()); + iSpkAngleRgt += spokeAngleStep; + iSpkAngleLft += spokeAngleStep; + } + + if ( version > 1 ) { + double maxHang = tgtPS->hubOverhangUS(); + double tmpHang = tgtPS->hubOverhangDS(); + if ( tmpHang > maxHang ) maxHang = tmpHang; + tgtPS->_envelHalfLength += maxHang + 2.0; + } else { + tgtPS->_envelHalfLength += Hub_overhang_Length; + } } else if (Hub_overhang_angle==90.0) { - HubRgtCornersZ.push_back(tgtPS->halfLength()); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); - - HubRgtCornersZ.push_back(tgtPS->halfLength()-Hub_thickness); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); - - HubRgtCornersZ.push_back(tgtPS->halfLength()-(Hub_thickness+0.000000001)); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); - - HubRgtCornersZ.push_back(tgtPS->halfLength()-Hub_hang_Length); - HubRgtCornersInnRadii.push_back(tgtPS->rOut()); - HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); - - HubLftCornersZ.push_back(-tgtPS->halfLength()); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); - - HubLftCornersZ.push_back(-tgtPS->halfLength()+Hub_thickness); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); - - HubLftCornersZ.push_back(-tgtPS->halfLength()+(Hub_thickness+0.000000001)); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); - - HubLftCornersZ.push_back(-tgtPS->halfLength()+Hub_hang_Length); - HubLftCornersInnRadii.push_back(tgtPS->rOut()); - HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); - - double SpkAnchrPosZ = HubRgtCornersZ.at(1); - double SpkAnchrPosRad = HubRgtCornersOutRadii.at(1)-spokeAnchordist; - double iSpkAngleRgt(0.0), iSpkAngleLft(-spokeSideDangle), tmpAngle(0.0); - for (int iSpk=0; iSpk::iterator pntPos; - pntPos = tgtPS->_anchoringPntsRgt.insert( - std::pair( - iSpkAngleRgt, - CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), - SpkAnchrPosRad*sin(tmpAngle), - SpkAnchrPosZ) - ) - ).first; - pntPos->second.transform(tgtPS->protonBeamRotation()); - tmpAngle = iSpkAngleLft * CLHEP::deg; - pntPos = tgtPS->_anchoringPntsLft.insert( - std::pair( - iSpkAngleLft, - CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), - SpkAnchrPosRad*sin(tmpAngle), - -SpkAnchrPosZ) - ) - ).first; - pntPos->second.transform(tgtPS->protonBeamRotation()); - iSpkAngleRgt += spokeAngleStep; - iSpkAngleLft += spokeAngleStep; - } + HubRgtCornersZ.push_back(tgtPS->halfLength()); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); + + HubRgtCornersZ.push_back(tgtPS->halfLength()-Hub_thickness); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); + + HubRgtCornersZ.push_back(tgtPS->halfLength()-(Hub_thickness+0.000000001)); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); + + HubRgtCornersZ.push_back(tgtPS->halfLength()-Hub_hang_Length); + HubRgtCornersInnRadii.push_back(tgtPS->rOut()); + HubRgtCornersOutRadii.push_back(HubRgtCornersInnRadii.back()+Hub_thickness); + + HubLftCornersZ.push_back(-tgtPS->halfLength()); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); + + HubLftCornersZ.push_back(-tgtPS->halfLength()+Hub_thickness); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness+Hub_overhang_Length); + + HubLftCornersZ.push_back(-tgtPS->halfLength()+(Hub_thickness+0.000000001)); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); + + HubLftCornersZ.push_back(-tgtPS->halfLength()+Hub_hang_Length); + HubLftCornersInnRadii.push_back(tgtPS->rOut()); + HubLftCornersOutRadii.push_back(HubLftCornersInnRadii.back()+Hub_thickness); + + double SpkAnchrPosZ = HubRgtCornersZ.at(1); + double SpkAnchrPosRad = HubRgtCornersOutRadii.at(1)-spokeAnchordist; + double iSpkAngleRgt(0.0), iSpkAngleLft(-spokeSideDangle), tmpAngle(0.0); + for (int iSpk=0; iSpk::iterator pntPos; + pntPos = tgtPS->_anchoringPntsRgt.insert( + std::pair( + iSpkAngleRgt, + CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), + SpkAnchrPosRad*sin(tmpAngle), + SpkAnchrPosZ) + ) + ).first; + pntPos->second.transform(tgtPS->protonBeamRotation()); + tmpAngle = iSpkAngleLft * CLHEP::deg; + pntPos = tgtPS->_anchoringPntsLft.insert( + std::pair( + iSpkAngleLft, + CLHEP::Hep3Vector(SpkAnchrPosRad*cos(tmpAngle), + SpkAnchrPosRad*sin(tmpAngle), + -SpkAnchrPosZ) + ) + ).first; + pntPos->second.transform(tgtPS->protonBeamRotation()); + iSpkAngleRgt += spokeAngleStep; + iSpkAngleLft += spokeAngleStep; + } } else { - throw cet::exception("GEOM") - << "Production Target, angle of the Hub overhang is not implemented \n"; + throw cet::exception("GEOM") + << "Production Target, angle of the Hub overhang is not implemented \n"; } tgtPS->_envelHalfLength *= cos(trgtMaxAngle); @@ -280,30 +280,30 @@ namespace mu2e { //const CLHEP::Hep3Vector& tgtPS_pos = tgtPS->position(); const CLHEP::Hep3Vector center(0,0,0); -// tgtPS->_pHubsRgtParams = std::unique_ptr -// (new Polycone(HubRgtCornersZ, -// HubRgtCornersInnRadii, -// HubRgtCornersOutRadii, -// center, -// c.getString("targetPS_Hub_materialName"))); -// -// tgtPS->_pHubsLftParams = std::unique_ptr -// (new Polycone(HubLftCornersZ, -// HubLftCornersInnRadii, -// HubLftCornersOutRadii, -// center, -// c.getString("targetPS_Hub_materialName"))); + // tgtPS->_pHubsRgtParams = std::unique_ptr + // (new Polycone(HubRgtCornersZ, + // HubRgtCornersInnRadii, + // HubRgtCornersOutRadii, + // center, + // c.getString("targetPS_Hub_materialName"))); + // + // tgtPS->_pHubsLftParams = std::unique_ptr + // (new Polycone(HubLftCornersZ, + // HubLftCornersInnRadii, + // HubLftCornersOutRadii, + // center, + // c.getString("targetPS_Hub_materialName"))); tgtPS->_pHubsRgtParams = new Polycone(HubRgtCornersZ, - HubRgtCornersInnRadii, - HubRgtCornersOutRadii, - center, - c.getString("targetPS_Hub_materialName")); + HubRgtCornersInnRadii, + HubRgtCornersOutRadii, + center, + c.getString("targetPS_Hub_materialName")); tgtPS->_pHubsLftParams = new Polycone(HubLftCornersZ, - HubLftCornersInnRadii, - HubLftCornersOutRadii, - center, - c.getString("targetPS_Hub_materialName")); + HubLftCornersInnRadii, + HubLftCornersOutRadii, + center, + c.getString("targetPS_Hub_materialName")); return tgtPS; @@ -339,9 +339,9 @@ namespace mu2e { c.getDouble("targetPS_rotY") * CLHEP::degree, c.getDouble("targetPS_rotZ") * CLHEP::degree, CLHEP::Hep3Vector(solenoidOffset, - 0, - c.getDouble("productionTarget.zNominal") - ) + 0, + c.getDouble("productionTarget.zNominal") + ) + c.getHep3Vector("productionTarget.offset"), c.getString("targetPS_targetCoreMaterial"), c.getString("targetPS_targetFinMaterial"), @@ -363,7 +363,7 @@ namespace mu2e { c.getDouble("targetPS_supportRingOuterRadius"), c.getDouble("targetPS_supportRingCutoutThickness"), c.getDouble("targetPS_supportRingCutoutLength") - )); + )); //check if we should configure supports //switch to using '.' instead of '_' to differentiate levels tgtPS->_supportsBuild = c.getBool("targetPS.supports.build", false); @@ -445,64 +445,64 @@ namespace mu2e { // Build parameter structs for Stickman constructor StickmanEnvelopeParams envelopeParams{ c.getDouble("targetPS_productionTargetMotherOuterRadius"), - c.getDouble("targetPS_productionTargetMotherHalfLength"), - c.getDouble("targetPS_rotX") * CLHEP::degree, - c.getDouble("targetPS_rotY") * CLHEP::degree, - c.getDouble("targetPS_rotZ") * CLHEP::degree, - c.getDouble("targetPS_halfStickmanLength"), - CLHEP::Hep3Vector(solenoidOffset, - 0, - c.getDouble("productionTarget.zNominal")) - + c.getHep3Vector("productionTarget.offset"), - c.getString("targetPS_targetVacuumMaterial") + c.getDouble("targetPS_productionTargetMotherHalfLength"), + c.getDouble("targetPS_rotX") * CLHEP::degree, + c.getDouble("targetPS_rotY") * CLHEP::degree, + c.getDouble("targetPS_rotZ") * CLHEP::degree, + c.getDouble("targetPS_halfStickmanLength"), + CLHEP::Hep3Vector(solenoidOffset, + 0, + c.getDouble("productionTarget.zNominal")) + + c.getHep3Vector("productionTarget.offset"), + c.getString("targetPS_targetVacuumMaterial") }; StickmanPlateParams plateParams{ nPlates, - plateMaterial, - plateROut, - nFins, - plateFinAngles, - c.getDouble("targetPS_plateFinOuterRadius"), - c.getDouble("targetPS_plateFinWidth"), - c.getDouble("targetPS_plateCenterToLugCenter"), - c.getDouble("targetPS_plateLugInnerRadius"), - c.getDouble("targetPS_plateLugOuterRadius"), - plateThickness, - plateLugThickness + plateMaterial, + plateROut, + nFins, + plateFinAngles, + c.getDouble("targetPS_plateFinOuterRadius"), + c.getDouble("targetPS_plateFinWidth"), + c.getDouble("targetPS_plateCenterToLugCenter"), + c.getDouble("targetPS_plateLugInnerRadius"), + c.getDouble("targetPS_plateLugOuterRadius"), + plateThickness, + plateLugThickness }; StickmanRodParams rodParams{ c.getString("targetPS_rodMaterial"), - c.getDouble("targetPS_rodRadius") + c.getDouble("targetPS_rodRadius") }; StickmanSpacerParams spacerParams{ c.getString("targetPS_spacerMaterial"), - c.getDouble("targetPS_spacerHalfLength"), - c.getDouble("targetPS_spacerOuterRadius"), - c.getDouble("targetPS_spacerInnerRadius") + c.getDouble("targetPS_spacerHalfLength"), + c.getDouble("targetPS_spacerOuterRadius"), + c.getDouble("targetPS_spacerInnerRadius") }; StickmanSupportRingParams supportRingParams{ c.getString("targetPS_supportRingMaterial"), - c.getDouble("targetPS_supportRingLength"), - c.getDouble("targetPS_supportRingInnerRadius"), - c.getDouble("targetPS_supportRingOuterRadius"), - c.getDouble("targetPS_supportRingLugOuterRadius"), - c.getDouble("targetPS_supportRingCutoutOffset") + c.getDouble("targetPS_supportRingLength"), + c.getDouble("targetPS_supportRingInnerRadius"), + c.getDouble("targetPS_supportRingOuterRadius"), + c.getDouble("targetPS_supportRingLugOuterRadius"), + c.getDouble("targetPS_supportRingCutoutOffset") }; std::unique_ptr tgtPS (new ProductionTarget( - c.getString("targetPS_model","NULL"), - c.getInt("targetPS_version"), - envelopeParams, - plateParams, - rodParams, - spacerParams, - supportRingParams - )); + c.getString("targetPS_model","NULL"), + c.getInt("targetPS_version"), + envelopeParams, + plateParams, + rodParams, + spacerParams, + supportRingParams + )); // Create configuration parameters struct StickmanConfigParams configParams; diff --git a/KinKalGeom/CMakeLists.txt b/KinKalGeom/CMakeLists.txt index c600163619..7bc59e97ac 100644 --- a/KinKalGeom/CMakeLists.txt +++ b/KinKalGeom/CMakeLists.txt @@ -2,7 +2,11 @@ cet_make_library( SOURCE src/CRV.cc src/KinKalGeom.cc + src/KKMaterial.cc + src/KKStrawMaterial.cc LIBRARIES PUBLIC + Offline::ConfigTools + KinKal::MatEnv KinKal::Geometry KinKal::Trajectory KinKal::General diff --git a/KinKalGeom/fcl/prolog.fcl b/KinKalGeom/fcl/prolog.fcl new file mode 100644 index 0000000000..ddf2db8258 --- /dev/null +++ b/KinKalGeom/fcl/prolog.fcl @@ -0,0 +1,18 @@ +BEGIN_PROLOG +KinKalGeom : { + KKMaterial : { + elements : "Offline/TrackerConditions/data/ElementsList.data" + isotopes : "Offline/TrackerConditions/data/IsotopesList.data" + materials : "Offline/TrackerConditions/data/MaterialsList.data" + strawGasMaterialName : "straw-gas" + strawWallMaterialName : "straw-wall" + strawWireMaterialName : "straw-wire" + IPAMaterialName : "HDPE" + STMaterialName : "Target" + IonizationEnergyLossMode : 1 # Moyal mean + SolidScatteringFraction : 0.999999 + GasScatteringFraction : 0.9999999 + ElectronBrehmsFraction : 0.04 + } +} +END_PROLOG diff --git a/KinKalGeom/inc/DetectorSolenoid.hh b/KinKalGeom/inc/DetectorSolenoid.hh index 7721710689..5fd98a55db 100644 --- a/KinKalGeom/inc/DetectorSolenoid.hh +++ b/KinKalGeom/inc/DetectorSolenoid.hh @@ -48,8 +48,8 @@ namespace mu2e { auto const& outerProtonAbsorberPtr() const { return opa_; } auto const& upstreamAbsorberPtr() const { return tsda_; } private: - CylPtr inner_; // inner cryostat cylinder - CylPtr outer_; // outer cryostat cylinder + CylPtr inner_; // inner cryostat cylinder boundary + CylPtr outer_; // outer cryostat cylinder boundary DiskPtr front_; // front (upstream) and back (downstream) of DS DiskPtr back_; CylPtr ipa_; // inner proton absorber diff --git a/Mu2eKinKal/inc/KKMaterial.hh b/KinKalGeom/inc/KKMaterial.hh similarity index 71% rename from Mu2eKinKal/inc/KKMaterial.hh rename to KinKalGeom/inc/KKMaterial.hh index 127d59642d..ef8a221ea8 100644 --- a/Mu2eKinKal/inc/KKMaterial.hh +++ b/KinKalGeom/inc/KKMaterial.hh @@ -1,5 +1,5 @@ -#ifndef Mu2eKinKal_KKMaterial_hh -#define Mu2eKinKal_KKMaterial_hh +#ifndef KinKalGeom_KKMaterial_hh +#define KinKalGeom_KKMaterial_hh // // build KinKal DetMaterial objects from art parameter configuration // @@ -8,14 +8,18 @@ #include "fhiclcpp/types/Tuple.h" // KinKal #include "KinKal/MatEnv/MatDBInfo.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" -#include "Offline/Mu2eKinKal/inc/KKFileFinder.hh" +#include "KinKal/MatEnv/FileFinderInterface.hh" +// KKGeom +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" +// mu2e +#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" +#include "Offline/Mu2eInterfaces/inc/Detector.hh" #include #include namespace mu2e { - class KKMaterial { + class KKMaterial : public MatEnv::FileFinderInterface, public Detector { public: using Name = fhicl::Name; using Comment = fhicl::Comment; @@ -24,7 +28,6 @@ namespace mu2e { fhicl::Atom isotopes { Name("isotopes"), Comment("Filename for istotopes information")}; fhicl::Atom elements { Name("elements"), Comment("Filename for elements information") }; fhicl::Atom materials { Name("materials"), Comment("Filename for materials information") }; - fhicl::Atom eloss { Name("ELossMode"), Comment("Energy Loss model (0=MPV, 1=Moyal"),MatEnv::DetMaterial::moyalmean }; fhicl::Atom strawGasMaterialName{ Name("strawGasMaterialName"), Comment("strawGasMaterialName") }; fhicl::Atom strawWallMaterialName{ Name("strawWallMaterialName"), Comment("strawWallMaterialName") }; fhicl::Atom strawWireMaterialName{ Name("strawWireMaterialName"), Comment("strawWireMaterialName") }; @@ -40,11 +43,22 @@ namespace mu2e { KKStrawMaterial const& strawMaterial() const; auto IPAMaterial() const { return matdbinfo_->findDetMaterial(ipamatname_); } auto STMaterial() const { return matdbinfo_->findDetMaterial(stmatname_); } + + // FileFinder interface + std::string matElmDictionaryFileName() const override; + std::string matIsoDictionaryFileName() const override; + std::string matMtrDictionaryFileName() const override; + std::string findFile( std::string const& path ) const override; private: - KKFileFinder filefinder_; // used to find material info + // material description files base names (not path) + std::string elementsBaseName_; + std::string isotopesBaseName_; + std::string materialsBaseName_; + mutable ConfigFileLookupPolicy policy_; + // specific material names std::string wallmatname_, gasmatname_, wirematname_,ipamatname_, stmatname_; mutable std::unique_ptr matdbinfo_; // material database - mutable std::unique_ptr smat_; // straw material + mutable std::unique_ptr smat_; // straw material; move to KinKalGeom }; } #endif diff --git a/Mu2eKinKal/inc/KKStrawMaterial.hh b/KinKalGeom/inc/KKStrawMaterial.hh similarity index 81% rename from Mu2eKinKal/inc/KKStrawMaterial.hh rename to KinKalGeom/inc/KKStrawMaterial.hh index ab76fe02b1..d29f153a12 100644 --- a/Mu2eKinKal/inc/KKStrawMaterial.hh +++ b/KinKalGeom/inc/KKStrawMaterial.hh @@ -6,7 +6,7 @@ // #include "KinKal/Trajectory/ClosestApproachData.hh" #include "Offline/TrackerGeom/inc/StrawProperties.hh" -#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" +#include namespace MatEnv { class MatDBInfo; @@ -31,14 +31,14 @@ namespace mu2e { const std::shared_ptr wiremat_); // construct using materials by name KKStrawMaterial(MatDBInfo const& matdbinfo,StrawProperties const& sprops, - const std::string& wallmat="straw-wall", const std::string& gasmat="straw-gas", const std::string& wiremat="straw-wire"); + const std::string& wallmat, const std::string& gasmat, const std::string& wiremat); // pathlength through straw components, given closest approach. Return the method used to compute the paths - PathCalc pathLengths(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, double& wallpath, double& gaspath, double& wirepath) const; + PathCalc pathLengths(ClosestApproachData const& cadata,double nsig, double& wallpath, double& gaspath, double& wirepath,int diag=0) const; PathCalc averagePathLengths(double& wallpath, double& gaspath, double& wirepath) const; // transit length given closest approach double transitLength(ClosestApproachData const& cadata) const; - // find the material crossings given doca and error on doca. Should allow for straw and wire to have different axes TODO - PathCalc findXings(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, std::vector& mxings) const; + // find the material crossings given doca and error on doca. The CA object should be WRT the straw axis, not the wire axis + PathCalc findXings(ClosestApproachData const& cadata, double nsig, std::vector& mxings, int diag=0) const; DetMaterial const& wallMaterial() const { return *wallmat_; } DetMaterial const& gasMaterial() const { return *gasmat_; } DetMaterial const& wireMaterial() const { return *wiremat_; } diff --git a/KinKalGeom/inc/KinKalGeom.hh b/KinKalGeom/inc/KinKalGeom.hh index 05770e6730..87c04e790e 100644 --- a/KinKalGeom/inc/KinKalGeom.hh +++ b/KinKalGeom/inc/KinKalGeom.hh @@ -13,20 +13,19 @@ #include "Offline/DataProducts/inc/SurfaceId.hh" #include "KinKal/Geometry/Surface.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" -#include "Offline/Mu2eInterfaces/inc/ProditionsEntity.hh" #include #include #include namespace mu2e { - class KinKalGeom : public Detector, public ProditionsEntity { + class KinKalGeom : public Detector { public: using SurfacePtr = std::shared_ptr; using SurfacePair =std::pair; using SurfacePairCollection = std::vector; using SurfacePairIter = std::multimap::const_iterator; using KKGMap = std::multimap; - // default constructor, now using GeometryService - KinKalGeom(); + // default constructor, now using GeometryService to create content + KinKalGeom(){} // accessor to the raw map auto const& map() const { return map_; } // find all surfaces that match an Id. Return vector can have >1 entry if the wildcard index (-1) is provided diff --git a/KinKalGeom/inc/ShellMaterial.hh b/KinKalGeom/inc/ShellMaterial.hh new file mode 100644 index 0000000000..91a8338949 --- /dev/null +++ b/KinKalGeom/inc/ShellMaterial.hh @@ -0,0 +1,21 @@ +// +// Class modeling a thin material as a surface, material, and thickness. Can be used to create a ShellXing in KKTrack +// Original author: David Brown (LBNL) 5/26 +// +#ifndef KinKalGeom_ShellMaterial_hh +#define KinKalGeom_ShellMaterial_hh +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "KinKal/Geometry/Surface.hh" +#include "KinKal/MatEnv/DetMaterial.hh" + +namespace mu2e { + class ShellMaterial { + using SurfacePtr = std::shared_ptr; + public: + private: + + SurfacePtr surf_; // surface for this materia + + + }; +} diff --git a/Mu2eKinKal/src/KKMaterial.cc b/KinKalGeom/src/KKMaterial.cc similarity index 67% rename from Mu2eKinKal/src/KKMaterial.cc rename to KinKalGeom/src/KKMaterial.cc index 126193fc9b..e80edefdb5 100644 --- a/Mu2eKinKal/src/KKMaterial.cc +++ b/KinKalGeom/src/KKMaterial.cc @@ -1,15 +1,16 @@ -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "KinKal/MatEnv/DetMaterial.hh" -#include "Offline/GeometryService/inc/GeometryService.hh" namespace mu2e { using MatDBInfo = MatEnv::MatDBInfo; using MatEnv::DetMaterial; KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : - filefinder_(matconfig.elements(),matconfig.isotopes(),matconfig.materials()), + elementsBaseName_(matconfig.elements()), + isotopesBaseName_(matconfig.isotopes()), + materialsBaseName_(matconfig.materials()), wallmatname_(matconfig.strawWallMaterialName()), gasmatname_(matconfig.strawGasMaterialName()), wirematname_(matconfig.strawWireMaterialName()), @@ -20,7 +21,7 @@ namespace mu2e { dmconf.scatterfrac_solid_ = matconfig.solidScatter(); dmconf.scatterfrac_gas_ = matconfig.gasScatter(); dmconf.ebrehmsfrac_ = matconfig.eBrehms(); - matdbinfo_ = std::make_unique(filefinder_,dmconf); + matdbinfo_ = std::make_unique(*this,dmconf); } KKStrawMaterial const& KKMaterial::strawMaterial() const { @@ -36,4 +37,10 @@ namespace mu2e { } return *smat_; } + + std::string KKMaterial::findFile( std::string const& basename ) const { return policy_( basename ); } + std::string KKMaterial::matElmDictionaryFileName() const { return findFile(elementsBaseName_); } + std::string KKMaterial::matIsoDictionaryFileName() const { return findFile(isotopesBaseName_); } + std::string KKMaterial::matMtrDictionaryFileName() const { return findFile(materialsBaseName_); } + } diff --git a/Mu2eKinKal/src/KKStrawMaterial.cc b/KinKalGeom/src/KKStrawMaterial.cc similarity index 88% rename from Mu2eKinKal/src/KKStrawMaterial.cc rename to KinKalGeom/src/KKStrawMaterial.cc index 4a7ee85d60..5c6e2f49cf 100644 --- a/Mu2eKinKal/src/KKStrawMaterial.cc +++ b/KinKalGeom/src/KKStrawMaterial.cc @@ -1,5 +1,4 @@ -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" -#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "KinKal/MatEnv/DetMaterial.hh" #include "KinKal/MatEnv/MatDBInfo.hh" #include "KinKal/Detector/MaterialXing.hh" @@ -35,11 +34,11 @@ namespace mu2e { return KKStrawMaterial::average; } - KKStrawMaterial::PathCalc KKStrawMaterial::pathLengths(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, - double& wallpath, double& gaspath, double& wirepath) const { + KKStrawMaterial::PathCalc KKStrawMaterial::pathLengths(ClosestApproachData const& cadata,double nsig, + double& wallpath, double& gaspath, double& wirepath,int diag) const { wallpath = gaspath = wirepath = 0.0; PathCalc retval = KKStrawMaterial::unknown; - double docarange = caconfig.nsig_*sqrt(std::max(0.0,cadata.docaVar())); + double docarange = nsig*sqrt(std::max(0.0,cadata.docaVar())); // if the doca range covers the straw, use the average double adoca = fabs(cadata.doca()); double mindoca = std::max(0.0,std::min(adoca,irad_)-docarange); @@ -58,10 +57,9 @@ namespace mu2e { if(awall>0.0)wallpath = awall/(omaxdoca-mindoca); retval = range; } - if(caconfig.diag_>0){ + if(diag>0){ std::cout << "KKStrawMaterial: DOCA " << fabs(cadata.doca()) << " DOCA Var " << cadata.docaVar() << " range [" << mindoca << "," << imaxdoca << "] wall path (mm)" << wallpath << " gas path " << gaspath << " retval " << retval << std::endl; } - // Model the wire as a diffuse gas, density constrained by DOCA TODO // 3D pathlength includes projection along straw double afac = angleFactor(cadata.dirDot()); @@ -79,15 +77,14 @@ namespace mu2e { return tlen; } - KKStrawMaterial::PathCalc KKStrawMaterial::findXings(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, - std::vector& mxings) const { + KKStrawMaterial::PathCalc KKStrawMaterial::findXings(ClosestApproachData const& cadata, double nsig, std::vector& mxings, int diag) const { mxings.clear(); double wallpath, gaspath, wirepath; - auto retval = pathLengths(cadata,caconfig,wallpath, gaspath, wirepath); + auto retval = pathLengths(cadata,nsig,wallpath, gaspath, wirepath,diag); if(wallpath > 0.0) mxings.push_back(MaterialXing(*wallmat_,wallpath)); if(gaspath > 0.0) mxings.push_back(MaterialXing(*gasmat_,gaspath)); if(wirepath > 0.0) mxings.push_back(MaterialXing(*wiremat_,wirepath)); - if(caconfig.diag_ > 1){ + if(diag > 1){ std::cout << "Ce wall dE/dx " << wallmat_->energyLoss(105,1.0,0.511) << std::endl; std::cout << "Ce gas dE/dx " << gasmat_->energyLoss(105,1.0,0.511) << std::endl; } diff --git a/KinKalGeom/src/KinKalGeom.cc b/KinKalGeom/src/KinKalGeom.cc index e6c21fd2cd..13bf8adcff 100644 --- a/KinKalGeom/src/KinKalGeom.cc +++ b/KinKalGeom/src/KinKalGeom.cc @@ -1,7 +1,6 @@ #include "Offline/KinKalGeom/inc/KinKalGeom.hh" #include "cetlib_except/exception.h" namespace mu2e { - KinKalGeom::KinKalGeom() : ProditionsEntity("KinKalGeom") {} void KinKalGeom::surfaces(SurfaceId const& id, SurfacePairCollection& surfs) const { // find all surfaces that match the input, including index wildcard diff --git a/KinKalGeom/src/SConscript b/KinKalGeom/src/SConscript index decefd3fbe..325b8e8f5c 100644 --- a/KinKalGeom/src/SConscript +++ b/KinKalGeom/src/SConscript @@ -14,6 +14,8 @@ rootlibs = env['ROOTLIBS'] helper = mu2e_helper(env) mainlib = helper.make_mainlib([ + 'mu2e_ConfigTools', + 'KinKal_MatEnv', 'KinKal_Geometry', 'KinKal_Trajectory', 'KinKal_General', @@ -22,6 +24,7 @@ mainlib = helper.make_mainlib([ 'art_Framework_Services_Registry', 'art_Utilities', 'canvas', + 'cetlib', 'cetlib_except', 'CLHEP', 'Core', diff --git a/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt b/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt index c764852f0a..9c2cb72103 100644 --- a/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt +++ b/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt @@ -17,8 +17,8 @@ vector productionTarget.offset = { 0., 0., 0.}; // Mother volume double targetPS_productionTargetMotherOuterRadius = 200.; //mm -double targetPS_productionTargetMotherHalfLength = 120.; //mm full target length 232.2 mm, - // allow the mother volume to be slightly larger +double targetPS_productionTargetMotherHalfLength = 120.; //mm full target length 232.2 mm, +// allow the mother volume to be slightly larger // Overall rotation double targetPS_rotX = 0.0; @@ -34,22 +34,22 @@ string targetPS_targetVacuumMaterial = "PSVacuum"; // Start with target plates. Chose to specify certain parameters for each plate, so can easily play with plate-by-plate dimensions later on int targetPS_numberOfPlates = 35; vector targetPS_plateMaterial = { // counting along mu2e +z direction, i.e. from the proton downstream side to the upstream side - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 5 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 10 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 15 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 20 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 25 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 30 - "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718" // 35 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 5 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 10 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 15 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 20 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 25 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 30 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718" // 35 }; // does not separate fin & core material here vector targetPS_rOut = { - 3.15, 3.15, 3.15, 3.15, 3.15, // 5 - 3.15, 3.15, 3.15, 3.15, 3.15, // 10 - 3.15, 3.15, 3.15, 3.15, 3.15, // 15 - 3.15, 3.15, 3.15, 3.15, 3.15, // 20 - 3.15, 3.15, 3.15, 3.15, 3.15, // 25 - 3.15, 3.15, 3.15, 3.15, 3.15, // 30 - 3.15, 3.15, 3.15, 3.15, 3.15 // 35 + 3.15, 3.15, 3.15, 3.15, 3.15, // 5 + 3.15, 3.15, 3.15, 3.15, 3.15, // 10 + 3.15, 3.15, 3.15, 3.15, 3.15, // 15 + 3.15, 3.15, 3.15, 3.15, 3.15, // 20 + 3.15, 3.15, 3.15, 3.15, 3.15, // 25 + 3.15, 3.15, 3.15, 3.15, 3.15, // 30 + 3.15, 3.15, 3.15, 3.15, 3.15 // 35 }; // in mm, core part radius, one value for each int targetPS_nStickmanFins = 3; // since the plates are installed on rods, the number of fins is always the same vector targetPS_plateFinAngles = {285.,165.,45.}; // degrees @@ -59,22 +59,22 @@ double targetPS_plateCenterToLugCenter = 19.6; // mm, center of p double targetPS_plateLugInnerRadius = 1.525; // mm, inner radius of the lug double targetPS_plateLugOuterRadius = 3.0; // mm, outer radius of the lug vector targetPS_plateThickness = { - 5.0, 5.0, 5.0, 5.0, 5.0, // 5 - 5.0, 5.0, 5.0, 5.0, 5.0, // 10 - 5.0, 5.0, 5.0, 5.0, 5.0, // 15 - 5.0, 5.0, 5.0, 5.0, 5.0, // 20 - 5.0, 5.0, 5.0, 5.0, 5.0, // 25 - 5.0, 5.0, 5.0, 5.0, 5.0, // 30 - 5.0, 5.0, 5.0, 5.0, 5.0 // 35 + 5.0, 5.0, 5.0, 5.0, 5.0, // 5 + 5.0, 5.0, 5.0, 5.0, 5.0, // 10 + 5.0, 5.0, 5.0, 5.0, 5.0, // 15 + 5.0, 5.0, 5.0, 5.0, 5.0, // 20 + 5.0, 5.0, 5.0, 5.0, 5.0, // 25 + 5.0, 5.0, 5.0, 5.0, 5.0, // 30 + 5.0, 5.0, 5.0, 5.0, 5.0 // 35 }; // in mm, one value for each plate, core and fin thickness vector targetPS_plateLugThickness = { - 6.0, 6.0, 6.0, 6.0, 6.0, // 5 - 6.0, 6.0, 6.0, 6.0, 6.0, // 10 - 6.0, 6.0, 6.0, 6.0, 6.0, // 15 - 6.0, 6.0, 6.0, 6.0, 6.0, // 20 - 6.0, 6.0, 6.0, 6.0, 6.0, // 25 - 6.0, 6.0, 6.0, 6.0, 6.0, // 30 - 6.0, 6.0, 6.0, 6.0, 6.0 // 35 + 6.0, 6.0, 6.0, 6.0, 6.0, // 5 + 6.0, 6.0, 6.0, 6.0, 6.0, // 10 + 6.0, 6.0, 6.0, 6.0, 6.0, // 15 + 6.0, 6.0, 6.0, 6.0, 6.0, // 20 + 6.0, 6.0, 6.0, 6.0, 6.0, // 25 + 6.0, 6.0, 6.0, 6.0, 6.0, // 30 + 6.0, 6.0, 6.0, 6.0, 6.0 // 35 }; // in mm, one value for each plate, note lugs are flush with the fins on the proton downstream side bool targetPS_addFilletToPlateCore = true; // whether to add fillet to the intersection between plate core and fins bool targetPS_addFilletToPlateLug = true; // whether to add fillet to the intersection between fins and lugs @@ -107,10 +107,10 @@ bool targetPS_addFilletToSupportRingLug = true; double targetPS_supportRingLugFilletRadius = 1.0; bool targetPS_addCutoutToSupportRing = true; // the six cylindrical cutouts at an angle double targetPS_supportRingCutoutOffset = 4.78; // mm, distance from the center of the cutout to the side of the end ring closer to the target center. - // This value is a must-have. If no cutout is added, this determines where the support wire is attached from the inner edge of the end ring. +// This value is a must-have. If no cutout is added, this determines where the support wire is attached from the inner edge of the end ring. int targetPS_nSupportRingCutouts = 3; vector targetPS_supportRingCutoutAngles = {15., 135., 255.}; // degrees - //{-45., 75., 195.}; these have spokes inside so not cutting out +//{-45., 75., 195.}; these have spokes inside so not cutting out double targetPS_supportRingCutoutInnerRadius = 2.06; // mm double targetPS_supportRingCutoutTilt = 20.0; //degrees, angle of the cutout w.r.t the radial direction @@ -129,9 +129,9 @@ vector targetPS.supports.features.rIns = {158.75, 158.75, 158.75}; //m vector targetPS.supports.rods.angles = {-45., 75., 195.0}; //degrees vector targetPS.supports.rods.halfLength = {70., 70., 70.}; //mm vector targetPS.supports.rods.offset = {-35., -14.80, 38.62}; //mm, changed according to F10122537_A_DWG1 - //rod center to the channel centers +//rod center to the channel centers vector targetPS.supports.rods.pinOffset = {3.81, 3.81, 3.81}; //mm, according to F10121763_A_DWG1, pinholes used to fix - //the rods are off-center, need to add pinOffset/cos(rotY) to offset above +//the rods are off-center, need to add pinOffset/cos(rotY) to offset above vector targetPS.supports.rods.radius = {9., 9., 9.}; //mm vector targetPS.supports.rods.radialOffset = {177.8, 177.8, 177.8}; //mm vector targetPS.supports.rods.wireOffset.downstream = {13., 13., 13.}; //mm, on the surface diff --git a/Mu2eG4/src/constructTargetPS.cc b/Mu2eG4/src/constructTargetPS.cc index 42e5625a3c..4771329060 100644 --- a/Mu2eG4/src/constructTargetPS.cc +++ b/Mu2eG4/src/constructTargetPS.cc @@ -59,20 +59,20 @@ using namespace std; namespace mu2e { void placeSubtractionVolumeBoxTubs( - const std::string& name //volume info for main object - ,const CLHEP::Hep3Vector& translation - ,VolumeInfo const& parent - ,G4Box* boxCutout //action volume object - ,G4Tubs* tubsCutout - ,const CLHEP::HepRotation* rotAngle - ,G4Material* material// finish nesting args - ,G4RotationMatrix* rotOverall - ,G4ThreeVector const& offset - ,int const copyNo - ,const G4Colour color - ,const std::string& lookupToken - ,const bool verbose = false - ) + const std::string& name //volume info for main object + ,const CLHEP::Hep3Vector& translation + ,VolumeInfo const& parent + ,G4Box* boxCutout //action volume object + ,G4Tubs* tubsCutout + ,const CLHEP::HepRotation* rotAngle + ,G4Material* material// finish nesting args + ,G4RotationMatrix* rotOverall + ,G4ThreeVector const& offset + ,int const copyNo + ,const G4Colour color + ,const std::string& lookupToken + ,const bool verbose = false + ) { VolumeInfo boxWithTubsCutoutInfo(name,translation,parent.centerInWorld); boxWithTubsCutoutInfo.solid = new G4SubtractionSolid(name,boxCutout,tubsCutout, rotOverall, offset); @@ -92,7 +92,7 @@ namespace mu2e { int targetVersion = _config.getInt("targetPS_version"); auto it = targetRegistry.find(targetVersion); - + if (it != targetRegistry.end()) { it->second(parent, _config); } else { @@ -109,7 +109,7 @@ namespace mu2e { AntiLeakRegistry & reg = _helper.antiLeakRegistry(); int verbosityLevel = _config.getInt("PS.verbosityLevel"); - + //G4Material* psVacVesselMaterial = findMaterialOrThrow(psVacVesselInnerParams.materialName()); verbosityLevel >0 && @@ -133,15 +133,15 @@ namespace mu2e { TubsParams prodTargetMotherParams( 0., _clamp_supWheel_rOut, envHalfLength); VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", - prodTargetMotherParams, - parent.logical->GetMaterial(), - 0, - tgt->position() - parent.centerInMu2e(), - parent, - 0, - G4Colour::Blue(), - "PS" - ); + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->position() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); if(verbosityLevel > 0) { G4cout << "inside clause 1 of " << __func__ << G4endl; G4cout << "target position and hall origin = " << tgt->position() << "\n" << _hallOriginInMu2e << " " << @@ -161,14 +161,14 @@ namespace mu2e { TubsParams suppWheelParams( _supWheel_trgtPS_rIn, _supWheel_trgtPS_rOut, _supWheel_trgtPS_halfLength); VolumeInfo suppWheelInfo = nestTubs( "ProductionTargetSupportWheel", - suppWheelParams, - _suppWheelMaterial, - 0, - _loclCenter, - prodTargetMotherInfo, - 0, - G4Colour::Gray() - ); + suppWheelParams, + _suppWheelMaterial, + 0, + _loclCenter, + prodTargetMotherInfo, + 0, + G4Colour::Gray() + ); double const _clamp_supWheel_rIn = _config.getDouble("clamp_supWheel_rIn"); G4Material* _clampSpWheelMaterial = findMaterialOrThrow(_config.getString("clamp_supWheel_material")); @@ -178,14 +178,14 @@ namespace mu2e { TubsParams clampSpWheelParams( _clamp_supWheel_rIn, _clamp_supWheel_rOut, _clamp_supWheel_halfLength); VolumeInfo clampSpWheelInfoR = nestTubs( "ClampSupportWheel_R", - clampSpWheelParams, - _clampSpWheelMaterial, - 0, - _clampPosR, - prodTargetMotherInfo, - 0, - G4Colour::Gray() - ); + clampSpWheelParams, + _clampSpWheelMaterial, + 0, + _clampPosR, + prodTargetMotherInfo, + 0, + G4Colour::Gray() + ); // G4ThreeVector _clampPosL(0.0,0.0,-(_supWheel_trgtPS_halfLength+_clamp_supWheel_halfLength)); // @@ -207,14 +207,14 @@ namespace mu2e { geomOptions->loadEntry( _config, "ProductionTarget", "targetPS" ); VolumeInfo prodTargetInfo = nestTubs( "ProductionTarget", - prodTargetParams, - prodTargetMaterial, - &tgt->productionTargetRotation(), - _loclCenter, - prodTargetMotherInfo, - 0, - G4Colour::Magenta() - ); + prodTargetParams, + prodTargetMaterial, + &tgt->productionTargetRotation(), + _loclCenter, + prodTargetMotherInfo, + 0, + G4Colour::Magenta() + ); if(verbosityLevel > 0) G4cout << "_loclCenter for Tier 1 target = " << _loclCenter << G4endl; // Now add fins for version 2 @@ -222,7 +222,7 @@ namespace mu2e { // Length of fins must be calculated: double finHalfLength = (tgt->halfLength() * 2.0 - tgt->hubDistUS() - - tgt->hubDistDS() - tgt->hubLenUS() - tgt->hubLenDS())/2.0; // on the inner side, + - tgt->hubDistDS() - tgt->hubLenUS() - tgt->hubLenDS())/2.0; // on the inner side, // adjacent to the main target body. // Use the steeper of the two hub angles to angle the ends of the fin @@ -232,9 +232,9 @@ namespace mu2e { double finHalfLengthOut = finHalfLength + tgt->finHeight()*std::cos(theAngle*CLHEP::degree); G4Trd * myTrd = new G4Trd("FinTrapezoid", - finHalfLength, finHalfLengthOut, - tgt->finThickness()/2.0, tgt->finThickness()/2.0, - tgt->finHeight()/2.0); + finHalfLength, finHalfLengthOut, + tgt->finThickness()/2.0, tgt->finThickness()/2.0, + tgt->finHeight()/2.0); // std::vector finDims = {tgt->finThickness()/2.0,tgt->finHeight()/2.0,finHalfLength}; double finZoff = (tgt->hubDistUS() - tgt->hubDistDS() + tgt->hubLenUS() - tgt->hubLenDS())/2.0; // z-offset for fin @@ -255,50 +255,50 @@ namespace mu2e { rotFin2->rotateX(M_PI/3.0); rotFin3->rotateX(5.0*M_PI/3.0); VolumeInfo fin1Vol("ProductionTargetFin1", - _loclCenter, - prodTargetMotherInfo.centerInWorld); + _loclCenter, + prodTargetMotherInfo.centerInWorld); VolumeInfo fin2Vol("ProductionTargetFin2", - _loclCenter, - prodTargetMotherInfo.centerInWorld); + _loclCenter, + prodTargetMotherInfo.centerInWorld); VolumeInfo fin3Vol("ProductionTargetFin3", - _loclCenter, - prodTargetMotherInfo.centerInWorld); + _loclCenter, + prodTargetMotherInfo.centerInWorld); fin1Vol.solid = myTrd; fin2Vol.solid = myTrd; fin3Vol.solid = myTrd; finishNesting( fin1Vol, - prodTargetMaterial, - rotFin1, - finOffset1, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); + prodTargetMaterial, + rotFin1, + finOffset1, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); finishNesting( fin2Vol, - prodTargetMaterial, - rotFin2, - finOffset2, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); + prodTargetMaterial, + rotFin2, + finOffset2, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); finishNesting( fin3Vol, - prodTargetMaterial, - rotFin3, - finOffset3, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); + prodTargetMaterial, + rotFin3, + finOffset3, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); } @@ -309,27 +309,27 @@ namespace mu2e { Polycone const & pHubRgtParams = *tgt->getHubsRgtPtr(); VolumeInfo prodTargetHubRgtInfo = nestPolycone("ProductionTargetHubRgt", - pHubRgtParams.getPolyconsParams(), - findMaterialOrThrow(pHubRgtParams.materialName()), - &tgt->productionTargetRotation(), - pHubRgtParams.originInMu2e(), - prodTargetMotherInfo, - 0, - G4Colour::Magenta(), - "ProductionTarget" - ); + pHubRgtParams.getPolyconsParams(), + findMaterialOrThrow(pHubRgtParams.materialName()), + &tgt->productionTargetRotation(), + pHubRgtParams.originInMu2e(), + prodTargetMotherInfo, + 0, + G4Colour::Magenta(), + "ProductionTarget" + ); Polycone const & pHubLftParams = *tgt->getHubsLftPtr(); VolumeInfo prodTargetHubLftInfo = nestPolycone("ProductionTargetHubLft", - pHubLftParams.getPolyconsParams(), - findMaterialOrThrow(pHubLftParams.materialName()), - &tgt->productionTargetRotation(), - pHubLftParams.originInMu2e(), - prodTargetMotherInfo, - 0, - G4Colour::Magenta(), - "ProductionTarget" - ); + pHubLftParams.getPolyconsParams(), + findMaterialOrThrow(pHubLftParams.materialName()), + &tgt->productionTargetRotation(), + pHubLftParams.originInMu2e(), + prodTargetMotherInfo, + 0, + G4Colour::Magenta(), + "ProductionTarget" + ); CLHEP::Hep3Vector zax(0,0,1); double spokeRad = 0.5*_config.getDouble("targetPS_Spoke_diameter"); @@ -377,19 +377,19 @@ namespace mu2e { iSpokeName<<"ProductionTargetSpokeRgt_"<supportRingMaterial()); TubsParams prodTargetMotherParams( 0. - ,tgt->productionTargetMotherOuterRadius() - ,tgt->productionTargetMotherHalfLength()); + ,tgt->productionTargetMotherOuterRadius() + ,tgt->productionTargetMotherHalfLength()); G4ThreeVector _loclCenter(0.0,0.0,0.0); G4RotationMatrix* targetRotation = reg.add(G4RotationMatrix(tgt->productionTargetRotation().inverse())); if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << "target rotation = " << *targetRotation << G4endl;} VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", - prodTargetMotherParams, - parent.logical->GetMaterial(), - 0, - tgt->haymanProdTargetPosition() - parent.centerInMu2e(), - parent, - 0, - G4Colour::Blue(), - "PS" - ); + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->haymanProdTargetPosition() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); if (verbosityLevel > 0){ G4cout << "target position and hall origin = " - << tgt->haymanProdTargetPosition() << "\n" << + << tgt->haymanProdTargetPosition() << "\n" << _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; } if (verbosityLevel > 2){ G4cout << __PRETTY_FUNCTION__ << "created prodTargetMotherInfo " - << tgt->productionTargetMotherOuterRadius() - << " " <productionTargetMotherHalfLength() << G4endl; + << tgt->productionTargetMotherOuterRadius() + << " " <productionTargetMotherHalfLength() << G4endl; } // // construct each section @@ -557,7 +557,7 @@ namespace mu2e { } G4Box* cutoutBox = reg.add(G4Box("cutoutBox",tgt->supportRingCutoutThickness()/2.+ magicOffset, - tgt->haymanFinThickness() + magicOffset,tgt->supportRingCutoutLength()/2.+ magicOffset)); + tgt->haymanFinThickness() + magicOffset,tgt->supportRingCutoutLength()/2.+ magicOffset)); // get real @@ -580,116 +580,116 @@ namespace mu2e { CLHEP::Hep3Vector currentStartingSegmentCenter = tgt->productionTargetRotation().inverse()*_startingSegmentCenter; VolumeInfo startingSegment = nestTubs(startName, - startingSegmentParams, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentStartingSegmentCenter, - prodTargetMotherInfo, - ithSection, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); + startingSegmentParams, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentStartingSegmentCenter, + prodTargetMotherInfo, + ithSection, + prodTargetVisible, + G4Colour::Yellow(), + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); //build fins around starting segment for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin) - { - const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); - if (verbosityLevel > 1) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - - double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; - G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); + { + const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); + if (verbosityLevel > 1) + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - // - // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. - // need extra length for visualization tool and overlap avoidance - G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); - ++finCopyNumber; + double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; + G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); - G4ThreeVector finTranslation = currentStartingSegmentCenter; - // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle - double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)) - ,distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTranslation = finTranslation + offsetRelativeToCore; - - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); - placeSubtractionVolumeBoxTubs(name - ,finTranslation - ,prodTargetMotherInfo - ,finBox - ,tubCutout - ,rotFinsG4.at(ithFin) - ,prodTargetFinMaterial - ,nullptr - ,downshift - ,finCopyNumber - ,G4Colour::Blue() - ,"PS" - ,verbosityLevel>1); - - // - // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. - // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. - const std::string nameTop = "ProductionTargetFinTopStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); - - double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; - double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; - if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); - G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); - ++finTopCopyNumber; - - if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " - << currentStartingSegmentCenter << " " << currentHalfStartingSegment - << " " << gapRadius << G4endl;} - G4ThreeVector finTopLocation = currentStartingSegmentCenter - + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfStartingSegment + gapRadius); - // - // for debugging - CLHEP::Hep3Vector offsetTop(0.,0.,0.); - // - // rotate the offset vector - G4ThreeVector finTopTranslation = finTopLocation + offsetTop; - double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; - CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*std::cos(currentFinAngles.at(ithFin)) - ,distanceTopFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTopTranslation = finTopTranslation + offsetTopRelativeToCore; - if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} - CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); - // - // these are oriented with gap radius along z. Recall still in mother volume so this axis is OK - G4RotationMatrix* tubTopRotation = reg.add(G4RotationMatrix(CLHEP::HepRotation::IDENTITY)); + // + // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. + // need extra length for visualization tool and overlap avoidance + G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); + ++finCopyNumber; - tubTopRotation->rotateY(M_PI/2.); + G4ThreeVector finTranslation = currentStartingSegmentCenter; + // + // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, + // both of the fin about the z axis and then the entire production target rotation angle + double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)) + ,distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTranslation = finTranslation + offsetRelativeToCore; + + CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); + placeSubtractionVolumeBoxTubs(name + ,finTranslation + ,prodTargetMotherInfo + ,finBox + ,tubCutout + ,rotFinsG4.at(ithFin) + ,prodTargetFinMaterial + ,nullptr + ,downshift + ,finCopyNumber + ,G4Colour::Blue() + ,"PS" + ,verbosityLevel>1); - VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); - finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); - finishNesting(finTopWithCutoutInfo - ,prodTargetFinMaterial - ,rotFinsG4.at(ithFin) - ,finTopTranslation - ,prodTargetMotherInfo.logical - ,finTopCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); - // - // this check also uses finWithCutout. Indirectly it's a way of forcing an overlap check or it won't compile, since finWithCutout is never used otherwise. - // if (doSurfaceCheck){checkForOverlaps(finTopWithCutout,_config,0);} + // + // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. + // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. + const std::string nameTop = "ProductionTargetFinTopStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); + + double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; + double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; + if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); + G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); + ++finTopCopyNumber; + + if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " + << currentStartingSegmentCenter << " " << currentHalfStartingSegment + << " " << gapRadius << G4endl;} + G4ThreeVector finTopLocation = currentStartingSegmentCenter + + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfStartingSegment + gapRadius); + // + // for debugging + CLHEP::Hep3Vector offsetTop(0.,0.,0.); + // + // rotate the offset vector + G4ThreeVector finTopTranslation = finTopLocation + offsetTop; + double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; + CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*std::cos(currentFinAngles.at(ithFin)) + ,distanceTopFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTopTranslation = finTopTranslation + offsetTopRelativeToCore; + if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} + CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); + // + // these are oriented with gap radius along z. Recall still in mother volume so this axis is OK + G4RotationMatrix* tubTopRotation = reg.add(G4RotationMatrix(CLHEP::HepRotation::IDENTITY)); + + tubTopRotation->rotateY(M_PI/2.); + + VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); + finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); + finishNesting(finTopWithCutoutInfo + ,prodTargetFinMaterial + ,rotFinsG4.at(ithFin) + ,finTopTranslation + ,prodTargetMotherInfo.logical + ,finTopCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); + // + // this check also uses finWithCutout. Indirectly it's a way of forcing an overlap check or it won't compile, since finWithCutout is never used otherwise. + // if (doSurfaceCheck){checkForOverlaps(finTopWithCutout,_config,0);} - } + } // // for very first starting section, we need to build the support ring @@ -699,16 +699,16 @@ namespace mu2e { if (ithSection==0){ std::string name = "ProductionTargetNegativeEndRing"; G4Tubs* supportRing = new G4Tubs(name, - tgt->supportRingInnerRadius(), - tgt->supportRingOuterRadius(), - tgt->supportRingLength()/2., - 0., - 2.*M_PI); + tgt->supportRingInnerRadius(), + tgt->supportRingOuterRadius(), + tgt->supportRingLength()/2., + 0., + 2.*M_PI); if (verbosityLevel > 0){ G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" - << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() - << ", " << tgt->supportRingLength()/2. << ")" << G4endl; + << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() + << ", " << tgt->supportRingLength()/2. << ")" << G4endl; } // move ring to end of target and then back by cutout size (signs for negative side) G4ThreeVector ringTranslation = currentStartingSegmentCenter @@ -736,9 +736,9 @@ namespace mu2e { ++boxCopyNumber; if (verbosityLevel > 6) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} if (verbosityLevel > 2){G4cout << " cutoutHeightAboveTarget = " << tgt->supportRingInnerRadius() << " " - << " " << tgt->supportRingCutoutThickness() << " " << cutoutHeightAboveTarget << G4endl;} + << " " << tgt->supportRingCutoutThickness() << " " << cutoutHeightAboveTarget << G4endl;} // // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, @@ -746,20 +746,20 @@ namespace mu2e { double distanceFromCenter = cutoutHeightAboveTarget; CLHEP::Hep3Vector boxShift(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)), - tgt->supportRingCutoutLength()/2. + magicOffset); + tgt->supportRingCutoutLength()/2. + magicOffset); if (verbosityLevel > 3){ G4cout << "production target rotation angle = " << tgt->productionTargetRotation().getTheta() << G4endl; G4cout << "fin dump: " << currentFinAngles.at(ithFin) << " " << distanceFromCenter << " " << offsetRelativeToCore << G4endl; G4cout << "ithfin; current fin angle; rotFin; finTranslation = " << ithFin << " " - << currentFinAngles.at(ithFin) << " " << *rotFinsG4.at(ithFin) << " " << ringTranslation << G4endl; + << currentFinAngles.at(ithFin) << " " << *rotFinsG4.at(ithFin) << " " << ringTranslation << G4endl; G4cout << "production target rotation = " << tgt->productionTargetRotation() << G4endl; } if (ithFin == 0){ - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); }else { - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); } } @@ -774,14 +774,14 @@ namespace mu2e { VolumeInfo ringWithCutoutNegative(name,ringTranslation,prodTargetMotherInfo.centerInWorld); ringWithCutoutNegative.solid = ringWithCutoutSolid.at(tgt->nHaymanFins() - 1); finishNesting(ringWithCutoutNegative - ,prodTargetSupportRingMaterial - ,rotRing - ,ringTranslation - ,prodTargetMotherInfo.logical - ,finCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + ,prodTargetSupportRingMaterial + ,rotRing + ,ringTranslation + ,prodTargetMotherInfo.logical + ,finCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); } @@ -807,26 +807,26 @@ namespace mu2e { _currentZ += currentHalfSegment + currentGap; if (verbosityLevel > 0){ G4cout << __PRETTY_FUNCTION__ << " ithSection , current Z for segment center is at " << ithSection << " " - << ithSegment << " " <<_currentZ << G4endl; + << ithSegment << " " <<_currentZ << G4endl; } _segmentCenter.setZ(_currentZ); G4ThreeVector currentSegmentCenter = tgt->productionTargetRotation().inverse()*_segmentCenter; //CLHEP::Hep3Vector currentSegmentCenter = tgt->productionTargetRotation()*_segmentCenter; // CLHEP::Hep3Vector currentSegmentCenter = _segmentCenter; VolumeInfo coreSegment = nestTubs(name, - segmentParams, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentSegmentCenter, - prodTargetMotherInfo, - coreCopyNumber, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); + segmentParams, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentSegmentCenter, + prodTargetMotherInfo, + coreCopyNumber, + prodTargetVisible, + G4Colour::Yellow(), + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); if (verbosityLevel > 2){G4cout << "segment center after volumeinfo = " << _segmentCenter << G4endl;} // @@ -834,52 +834,52 @@ namespace mu2e { for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin) // for (int ithFin = 0; ithFin < 1; ++ithFin) - { - const std::string name = "ProductionTargetFinSection_" + std::to_string(ithSection+1) + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); - if (verbosityLevel > 1) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - // divide by 2 to make box half-height since I gave it the "radius" to start with; arbitrary convention. - double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; - if (verbosityLevel > 2){G4cout << "finHeightAboveTarget = " << tgt->finOuterRadius() << " " - << tgt->rOut() << " " << finHalfHeightAboveTarget << G4endl;} - G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfSegment); - // - // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. - G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfSegment + magicOffset,dSphi,dPphi);// need extra length for visualization tool - // - // I now have two G4Solids, a box and a cutout. Combine them to make a logical volume. the box is - // oriented so that it is "z" thick and "radius" high. the tubs is made to be the same, so no rotation - // and also no translation needed - - ++finCopyNumber; - - // - //for debugging - CLHEP::Hep3Vector offset(0.,0.,0.); - G4ThreeVector finTranslation = currentSegmentCenter; - - double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTranslation = finTranslation + offsetRelativeToCore; - - - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); - placeSubtractionVolumeBoxTubs(name - ,finTranslation - ,prodTargetMotherInfo - ,finBox - ,tubCutout - ,rotFinsG4.at(ithFin) - ,prodTargetFinMaterial - ,nullptr - ,downshift - ,finCopyNumber - ,G4Colour::Blue() - ,"PS" - ,verbosityLevel>1); + { + const std::string name = "ProductionTargetFinSection_" + std::to_string(ithSection+1) + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); + if (verbosityLevel > 1) + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} + // divide by 2 to make box half-height since I gave it the "radius" to start with; arbitrary convention. + double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; + if (verbosityLevel > 2){G4cout << "finHeightAboveTarget = " << tgt->finOuterRadius() << " " + << tgt->rOut() << " " << finHalfHeightAboveTarget << G4endl;} + G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfSegment); + // + // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. + G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfSegment + magicOffset,dSphi,dPphi);// need extra length for visualization tool + // + // I now have two G4Solids, a box and a cutout. Combine them to make a logical volume. the box is + // oriented so that it is "z" thick and "radius" high. the tubs is made to be the same, so no rotation + // and also no translation needed + + ++finCopyNumber; + + // + //for debugging + CLHEP::Hep3Vector offset(0.,0.,0.); + G4ThreeVector finTranslation = currentSegmentCenter; + + double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTranslation = finTranslation + offsetRelativeToCore; + + + CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); + placeSubtractionVolumeBoxTubs(name + ,finTranslation + ,prodTargetMotherInfo + ,finBox + ,tubCutout + ,rotFinsG4.at(ithFin) + ,prodTargetFinMaterial + ,nullptr + ,downshift + ,finCopyNumber + ,G4Colour::Blue() + ,"PS" + ,verbosityLevel>1); // // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. @@ -895,10 +895,10 @@ namespace mu2e { ++finTopCopyNumber; if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " - << currentStartingSegmentCenter << " " << currentHalfStartingSegment << " " - << gapRadius << G4endl;} + << currentStartingSegmentCenter << " " << currentHalfStartingSegment << " " + << gapRadius << G4endl;} G4ThreeVector finTopBoxLocation = currentSegmentCenter - + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfSegment + gapRadius); + + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfSegment + gapRadius); // // for debugging CLHEP::Hep3Vector offsetTop(0.,0.,0.); @@ -926,16 +926,16 @@ namespace mu2e { VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); finishNesting(finTopWithCutoutInfo - ,prodTargetFinMaterial - ,rotFinsG4.at(ithFin) - ,finTopTranslation - ,prodTargetMotherInfo.logical - ,finTopCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + ,prodTargetFinMaterial + ,rotFinsG4.at(ithFin) + ,finTopTranslation + ,prodTargetMotherInfo.logical + ,finTopCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); - } + } _currentZ += currentHalfSegment; // // if this is the last segment, add another gap before next section starts! @@ -946,7 +946,7 @@ namespace mu2e { } if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ending at z=" << _currentZ << G4endl;} } - ++coreCopyNumber; + ++coreCopyNumber; // // if this is the last section, add on one more beginning block. decided not to extend vectors and have zero segments, just ugly if (ithSection == (numberOfSections -1) ){ @@ -961,19 +961,19 @@ namespace mu2e { CLHEP::Hep3Vector currentStartingSegmentCenter = tgt->productionTargetRotation().inverse()*_startingSegmentCenter; double currentHalfStartingSegment = tgt->startingSectionThickness().at(ithSection)/2.; startingSegment = nestTubs(startName, - startingSegmentParamsEnd, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentStartingSegmentCenter, - prodTargetMotherInfo, - ithSection+1, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); + startingSegmentParamsEnd, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentStartingSegmentCenter, + prodTargetMotherInfo, + ithSection+1, + prodTargetVisible, + G4Colour::Yellow(), + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); //build fins around starting segment for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ @@ -982,7 +982,7 @@ namespace mu2e { const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1+1) + "_Fin_" + std::to_string(ithFin); if (verbosityLevel > 1) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); @@ -1006,8 +1006,8 @@ namespace mu2e { // // for debugging CLHEP::Hep3Vector offset(0.,0.,0.); - // - // rotate the offset vector + // + // rotate the offset vector G4ThreeVector finTranslation = currentStartingSegmentCenter; // // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, @@ -1024,14 +1024,14 @@ namespace mu2e { VolumeInfo finWithCutoutInfo(name,finTranslation,prodTargetMotherInfo.centerInWorld); finWithCutoutInfo.solid = new G4SubtractionSolid(name,finBox,tubCutout,nullptr,downshift); finishNesting(finWithCutoutInfo - ,prodTargetFinMaterial - ,rotFinsG4.at(ithFin) - ,finTranslation - ,prodTargetMotherInfo.logical - ,finCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + ,prodTargetFinMaterial + ,rotFinsG4.at(ithFin) + ,finTranslation + ,prodTargetMotherInfo.logical + ,finCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); } // @@ -1044,15 +1044,15 @@ namespace mu2e { /***********and now the final end ring ***********/ std::string name = "ProductionTargetPositiveEndRing"; G4Tubs* supportRing = new G4Tubs(name, - tgt->supportRingInnerRadius(), - tgt->supportRingOuterRadius(), - tgt->supportRingLength()/2., - 0., - 2.*M_PI); + tgt->supportRingInnerRadius(), + tgt->supportRingOuterRadius(), + tgt->supportRingLength()/2., + 0., + 2.*M_PI); if (verbosityLevel > 0){ G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" - << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() - << ", " << tgt->supportRingLength()/2. << ")" << G4endl; + << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() + << ", " << tgt->supportRingLength()/2. << ")" << G4endl; } // move ring to end of target and then back by cutout size (signs for positive side) G4ThreeVector ringTranslation = currentStartingSegmentCenter @@ -1079,7 +1079,7 @@ namespace mu2e { CLHEP::Hep3Vector boxShift(distanceFromCenter*std::cos(currentFinAngle),distanceFromCenter*std::sin(currentFinAngle),0.); CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngle),distanceFromCenter*std::sin(currentFinAngle), - -tgt->supportRingCutoutLength()/2. - magicOffset);// note flipped sign from other end of target + -tgt->supportRingCutoutLength()/2. - magicOffset);// note flipped sign from other end of target CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,0.,0.);//for debugging offsetRelativeToCore = offsetRelativeToCore + downshift; @@ -1100,14 +1100,14 @@ namespace mu2e { VolumeInfo ringWithCutoutPositive(name,ringTranslation,prodTargetMotherInfo.centerInWorld); ringWithCutoutPositive.solid = ringWithCutoutSolid.at(tgt->nHaymanFins() - 1); // what does this = mean? finishNesting(ringWithCutoutPositive - ,prodTargetSupportRingMaterial - ,rotRing - ,ringTranslation - ,prodTargetMotherInfo.logical - ,finCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + ,prodTargetSupportRingMaterial + ,rotRing + ,ringTranslation + ,prodTargetMotherInfo.logical + ,finCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); } //end if(ithSection == (numberOfSections -1) ) } //end for(int ithSection...) @@ -1120,16 +1120,16 @@ namespace mu2e { //create the volume info for the support wheel+rods VolumeInfo suppWheelInfo( "ProductionTargetSupportWheel", localWheelCenter, prodTargetMotherInfo.centerInMu2e()); suppWheelInfo.solid = new G4Tubs("ProductionTargetSupportWheel_wheel", suppWheelParams[0], suppWheelParams[1], - suppWheelParams[2], 0., CLHEP::twopi); - // suppWheelParams, - // suppWheelMaterial, - // 0, - // localWheelCenter, - // prodTargetMotherInfo, - // 0, - // G4Colour::Gray(), - // "PS" - // ); + suppWheelParams[2], 0., CLHEP::twopi); + // suppWheelParams, + // suppWheelMaterial, + // 0, + // localWheelCenter, + // prodTargetMotherInfo, + // 0, + // G4Colour::Gray(), + // "PS" + // ); // add spokes // @@ -1184,14 +1184,14 @@ namespace mu2e { CLHEP::Hep3Vector featureCenter(localWheelCenter); //center is wheel center double featureParams[] = {featureRIn, featureROut, tgt->supportWheelHL(), featureAngle - featureArc/2. /*phi0*/, featureArc /*dphi*/}; G4Tubs* featureTubs = new G4Tubs("ProductionTargetSupportFeature_" +std::to_string(ispoke), - featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); + featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelFeature_union_"+std::to_string(ispoke), - suppWheelInfo.solid, featureTubs, 0, featureCenter); + suppWheelInfo.solid, featureTubs, 0, featureCenter); //add the support rod to the wheel G4Tubs* rodTubs = new G4Tubs("ProductionTargetSupportRod_" +std::to_string(ispoke), - 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); + 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelRod_union_"+std::to_string(ispoke), - suppWheelInfo.solid, rodTubs, rodRot, rodCenter); + suppWheelInfo.solid, rodTubs, rodRot, rodCenter); } const int side = (1-2*istream); //+1 or -1 //info about wire connection @@ -1201,8 +1201,8 @@ namespace mu2e { wheelPos += side*supportWheelRodHL[ispoke]*rodAxis; //translate from rod center to edge CLHEP::Hep3Vector rodCenterToWire(std::cos(wheelAngle)*std::cos(targetAngle), - std::sin(wheelAngle)*std::cos(targetAngle), - -std::cos(wheelAngle)*std::sin(targetAngle)); + std::sin(wheelAngle)*std::cos(targetAngle), + -std::cos(wheelAngle)*std::sin(targetAngle)); wheelPos -= supportWheelRodRadius[ispoke]*rodCenterToWire; double zWireRodOffset = (istream == 0) ? supportWheelRodWireOffsetD[ispoke] : supportWheelRodWireOffsetU[ispoke]; wheelPos -= side*zWireRodOffset*rodAxis; @@ -1216,9 +1216,9 @@ namespace mu2e { CLHEP::Hep3Vector zax(0.,0.,1.); if(verbosityLevel > 0) std::cout << "istream " << istream << " ispoke " << ispoke << std::endl - << "target pos " << targetPos << "\nwheel pos " << wheelPos << std::endl - << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << std::endl - << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << std::endl; + << "target pos " << targetPos << "\nwheel pos " << wheelPos << std::endl + << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << std::endl + << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << std::endl; //to remove overlaps where the wire connects, need angle of wire and surface connecting to //remove overlap at target double wireTargetAngle = targetAxis.angle(-1.*spokeAxis); @@ -1226,7 +1226,7 @@ namespace mu2e { targetPos += (deltaLength+0.1)*spokeAxis; //subtract off the length if(verbosityLevel > 0) std::cout << "wire target angle " << wireTargetAngle << " delta L " << deltaLength - << " target pos " << targetPos < 0) std::cout << "wire rod angle " << wireRodAngle << " delta L " << deltaLength - << " wheel pos " << wheelPos <stickmanSupportRingMaterial()); TubsParams prodTargetMotherParams( 0. - ,tgt->productionTargetMotherOuterRadius() - ,tgt->productionTargetMotherHalfLength()); + ,tgt->productionTargetMotherOuterRadius() + ,tgt->productionTargetMotherHalfLength()); G4ThreeVector _localCenter(0.0,0.0,0.0); G4RotationMatrix* targetRotation = reg.add(G4RotationMatrix(tgt->productionTargetRotation().inverse())); if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << "target rotation = " << *targetRotation << G4endl;} VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", - prodTargetMotherParams, - parent.logical->GetMaterial(), - 0, - tgt->stickmanProdTargetPosition() - parent.centerInMu2e(), - parent, - 0, - G4Colour::Blue(), - "PS" - ); + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->stickmanProdTargetPosition() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); if (verbosityLevel > 0){ G4cout << "target position and hall origin = " - << tgt->stickmanProdTargetPosition() << "\n" - << _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; + << tgt->stickmanProdTargetPosition() << "\n" + << _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; } if (verbosityLevel > 2){ G4cout << __PRETTY_FUNCTION__ << "created prodTargetMotherInfo " - << tgt->productionTargetMotherOuterRadius() - << " " <productionTargetMotherHalfLength() << G4endl; + << tgt->productionTargetMotherOuterRadius() + << " " <productionTargetMotherHalfLength() << G4endl; } // tiny offset to avoid precision issues @@ -1378,7 +1378,7 @@ namespace mu2e { if (cached != stickmanPlateSolidCache.end()) { if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " reusing cached plate solid for plate " - << ithPlate << G4endl; + << ithPlate << G4endl; } return cached->second; } @@ -1391,41 +1391,41 @@ namespace mu2e { const double lugZShift = lugHalfThickness - plateHalfThickness; if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " creating new plate solid for plate " - << ithPlate << " with cache key " << cacheKey << G4endl; + << ithPlate << " with cache key " << cacheKey << G4endl; } // solid for the plate core G4VSolid* plateSolid = reg.add(new G4Tubs(cacheKey + "_Core", // - 0., // inner radius - tgt->plateROut(ithPlate), // outer radius - plateHalfThickness, // z half length - 0., // starting angle - CLHEP::twopi)); // segment angle is full circle + 0., // inner radius + tgt->plateROut(ithPlate), // outer radius + plateHalfThickness, // z half length + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " core solid parameters (rOut, halfLength) = (" - << tgt->plateROut(ithPlate) << ", " << plateHalfThickness << ")" << G4endl; + << tgt->plateROut(ithPlate) << ", " << plateHalfThickness << ")" << G4endl; } // solid for the fin G4Box* finBox = reg.add(new G4Box(cacheKey + "_FinBox", // - finHalfLength, // x half length - finHalfWidth, // y half length - plateHalfThickness)); // z half length + finHalfLength, // x half length + finHalfWidth, // y half length + plateHalfThickness)); // z half length if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " fin solid parameters (halfLength, halfWidth, halfThickness) = (" - << finHalfLength << ", " << finHalfWidth << ", " << plateHalfThickness << ")" << G4endl; + << finHalfLength << ", " << finHalfWidth << ", " << plateHalfThickness << ")" << G4endl; } // solid for the lug G4Tubs* lugSolid = reg.add(new G4Tubs(cacheKey + "_Lug", // - tgt->plateLugInnerRadius(), // inner radius - tgt->plateLugOuterRadius(), // outer radius - lugHalfThickness, // z half length - 0., // starting angle - CLHEP::twopi)); // segment angle is full circle + tgt->plateLugInnerRadius(), // inner radius + tgt->plateLugOuterRadius(), // outer radius + lugHalfThickness, // z half length + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " lug solid parameters (rIn, rOut, halfThickness) = (" - << tgt->plateLugInnerRadius() << ", " << tgt->plateLugOuterRadius() << ", " << lugHalfThickness << ")" << G4endl; + << tgt->plateLugInnerRadius() << ", " << tgt->plateLugOuterRadius() << ", " << lugHalfThickness << ")" << G4endl; } // add fillets to the plate core @@ -1446,7 +1446,7 @@ namespace mu2e { if (tgt->addFilletToPlateCore()) { // get the coordinates of the fillet union vertices for the fin at 0 degree const double coreFilletX = std::sqrt(std::pow(tgt->plateROut(ithPlate) + plateFilletRadius, 2) - - std::pow(finHalfWidth + plateFilletRadius, 2)); + - std::pow(finHalfWidth + plateFilletRadius, 2)); std::vector coreFilletVertices = { G4TwoVector(0., 0.), G4TwoVector(coreFilletX, finHalfWidth + plateFilletRadius), @@ -1454,46 +1454,46 @@ namespace mu2e { }; // create the extrusion solid for the fillet G4ExtrudedSolid* coreFilletExtrusion = reg.add(new G4ExtrudedSolid( - cacheKey + "_CoreFilletExtrusion", // - coreFilletVertices, // vertices - plateHalfThickness, // half length in z - G4TwoVector(0.,0.), // shift of polygon center at -half length - 1.0, // scale of polygon at -half length - G4TwoVector(0.,0.), // shift of polygon center at +half length - 1.0)); // scale of polygon at +half length + cacheKey + "_CoreFilletExtrusion", // + coreFilletVertices, // vertices + plateHalfThickness, // half length in z + G4TwoVector(0.,0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0.,0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " core fillet extrusion vertices = (" - << coreFilletVertices.at(0) << ", " - << coreFilletVertices.at(1) << ", " - << coreFilletVertices.at(2) << ")" << G4endl; + << coreFilletVertices.at(0) << ", " + << coreFilletVertices.at(1) << ", " + << coreFilletVertices.at(2) << ")" << G4endl; G4cout << __PRETTY_FUNCTION__ << " core fillet extrusion half length = " << plateHalfThickness << G4endl; } // create the cutout cylinder for the fillet G4Tubs* coreFilletCutout = reg.add(new G4Tubs( - cacheKey + "_CoreFilletCutout", // - 0., // inner radius - plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) - plateHalfThickness + stickmanMagicOffset, // half length in z - 0., // starting angle - CLHEP::twopi)); // segment angle is full circle + cacheKey + "_CoreFilletCutout", // + 0., // inner radius + plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle // make 2 cutouts for the fillet, one at the upper edge of the fin and one at the lower edge of the fin G4VSolid* coreFilletUpper = reg.add(new G4SubtractionSolid( - cacheKey + "_CoreFilletUpper", // - coreFilletExtrusion, // base solid - coreFilletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(coreFilletX, - finHalfWidth + plateFilletRadius, - 0.))); // translation of cutout solid + cacheKey + "_CoreFilletUpper", // + coreFilletExtrusion, // base solid + coreFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(coreFilletX, + finHalfWidth + plateFilletRadius, + 0.))); // translation of cutout solid plateCoreFillet = reg.add(new G4SubtractionSolid( - cacheKey + "_CoreFilletLower", // - coreFilletUpper, // base solid - coreFilletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(coreFilletX, - -finHalfWidth - plateFilletRadius, - 0.))); // translation of cutout solid + cacheKey + "_CoreFilletLower", // + coreFilletUpper, // base solid + coreFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(coreFilletX, + -finHalfWidth - plateFilletRadius, + 0.))); // translation of cutout solid if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " core fillet cutout radius = " << plateFilletRadius + stickmanMagicOffset << G4endl; G4cout << __PRETTY_FUNCTION__ << " core fillet upper cutout center = (" << coreFilletX << ", " << finHalfWidth + plateFilletRadius << ", 0)" << G4endl; @@ -1521,7 +1521,7 @@ namespace mu2e { // get the coordinates of the fillet union vertices for the fin at 0 degree const double lugFilletX = tgt->plateCenterToLugCenter() - std::sqrt(std::pow(tgt->plateLugOuterRadius() + plateFilletRadius, 2) - - std::pow(finHalfWidth + plateFilletRadius, 2)); + - std::pow(finHalfWidth + plateFilletRadius, 2)); std::vector lugFilletVertices = { G4TwoVector(tgt->plateCenterToLugCenter(), 0.), G4TwoVector(lugFilletX, finHalfWidth + plateFilletRadius), @@ -1529,62 +1529,62 @@ namespace mu2e { }; // create the extrusion solid for the fillet G4ExtrudedSolid* lugFilletExtrusion = reg.add(new G4ExtrudedSolid( - cacheKey + "_LugFilletExtrusion", // - lugFilletVertices, // vertices - plateHalfThickness, // half length in z - G4TwoVector(0.,0.), // shift of polygon center at -half length - 1.0, // scale of polygon at -half length - G4TwoVector(0.,0.), // shift of polygon center at +half length - 1.0)); // scale of polygon at +half length + cacheKey + "_LugFilletExtrusion", // + lugFilletVertices, // vertices + plateHalfThickness, // half length in z + G4TwoVector(0.,0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0.,0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " lug fillet extrusion vertices = (" - << lugFilletVertices.at(0) << ", " - << lugFilletVertices.at(1) << ", " - << lugFilletVertices.at(2) << ")" << G4endl; + << lugFilletVertices.at(0) << ", " + << lugFilletVertices.at(1) << ", " + << lugFilletVertices.at(2) << ")" << G4endl; G4cout << __PRETTY_FUNCTION__ << " lug fillet extrusion half length = " << plateHalfThickness << G4endl; } // create the cutout cylinder for the fillet G4Tubs* lugFilletCutout = reg.add(new G4Tubs( - cacheKey + "_LugFilletCutout", // - 0., // inner radius - plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) - plateHalfThickness + stickmanMagicOffset, // half length in z - 0., // starting angle - CLHEP::twopi)); // segment angle is full circle + cacheKey + "_LugFilletCutout", // + 0., // inner radius + plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle // lug center needs to be cut out too. Create a second cutout cylinder for the lug center with radius of // plateLugInnerRadius and height of plateThickness(ithPlate) centered at (plateCenterToLugCenter, 0) G4Tubs* lugCenterCutout = reg.add(new G4Tubs( - cacheKey + "_LugCenterCutout", // - 0., // inner radius - tgt->plateLugInnerRadius() + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) - plateHalfThickness + stickmanMagicOffset, // half length in z - 0., // starting angle - CLHEP::twopi)); // segment angle is full circle + cacheKey + "_LugCenterCutout", // + 0., // inner radius + tgt->plateLugInnerRadius() + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle // cut out the lug center from the extrusion first G4VSolid* lugFilletWithCenterCutout = reg.add(new G4SubtractionSolid( - cacheKey + "_LugFilletWithCenterCutout", // - lugFilletExtrusion, // base solid - lugCenterCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(tgt->plateCenterToLugCenter(), 0., 0.))); // translation of cutout solid + cacheKey + "_LugFilletWithCenterCutout", // + lugFilletExtrusion, // base solid + lugCenterCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(tgt->plateCenterToLugCenter(), 0., 0.))); // translation of cutout solid // make 2 cutouts for the fillet, one at the upper edge of the fin and one at the lower edge of the fin G4VSolid* lugFilletUpper = reg.add(new G4SubtractionSolid( - cacheKey + "_LugFilletUpper", // - lugFilletWithCenterCutout, // base solid - lugFilletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(lugFilletX, - finHalfWidth + plateFilletRadius, - 0.))); // translation of cutout solid + cacheKey + "_LugFilletUpper", // + lugFilletWithCenterCutout, // base solid + lugFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(lugFilletX, + finHalfWidth + plateFilletRadius, + 0.))); // translation of cutout solid plateLugFillet = reg.add(new G4SubtractionSolid( - cacheKey + "_LugFilletLower", // - lugFilletUpper, // base solid - lugFilletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(lugFilletX, - -finHalfWidth - plateFilletRadius, - 0.))); // translation of cutout solid + cacheKey + "_LugFilletLower", // + lugFilletUpper, // base solid + lugFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(lugFilletX, + -finHalfWidth - plateFilletRadius, + 0.))); // translation of cutout solid if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " lug fillet cutout radius = " << plateFilletRadius + stickmanMagicOffset << G4endl; G4cout << __PRETTY_FUNCTION__ << " lug center cutout radius = " << tgt->plateLugInnerRadius() + stickmanMagicOffset << G4endl; @@ -1599,19 +1599,19 @@ namespace mu2e { const double currentFinAngle = tgt->plateFinAngle(ithFin); // fin and lug shifts of the current fin const CLHEP::Hep3Vector finShift(finHalfLength * std::cos(currentFinAngle), - finHalfLength * std::sin(currentFinAngle), - 0.); + finHalfLength * std::sin(currentFinAngle), + 0.); const CLHEP::Hep3Vector lugShift(tgt->plateCenterToLugCenter() * std::cos(currentFinAngle), - tgt->plateCenterToLugCenter() * std::sin(currentFinAngle), - lugZShift); + tgt->plateCenterToLugCenter() * std::sin(currentFinAngle), + lugZShift); // note the lug is shifted in z to be flush with the edge of the plate core // add fin plateSolid = reg.add(new G4UnionSolid(cacheKey + "_FinUnion_" + std::to_string(ithFin), - plateSolid, // base solid - finBox, // union solid - stickmanFinRotations.at(ithFin), // rotation of union solid - finShift)); // translation of union solid + plateSolid, // base solid + finBox, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + finShift)); // translation of union solid if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " fin " << ithFin << " added with angle of " << tgt->plateFinAngle(ithFin)/CLHEP::degree << " deg" << G4endl; } @@ -1619,18 +1619,18 @@ namespace mu2e { // add core fillets if applicable if (plateCoreFillet != nullptr) { plateSolid = reg.add(new G4UnionSolid(cacheKey + "_CoreFilletUnion_" + std::to_string(ithFin), - plateSolid, // base solid - plateCoreFillet, // union solid - stickmanFinRotations.at(ithFin), // rotation of union solid - CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid + plateSolid, // base solid + plateCoreFillet, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid } // add lug plateSolid = reg.add(new G4UnionSolid(cacheKey + "_LugUnion_" + std::to_string(ithFin), - plateSolid, // base solid - lugSolid, // union solid - nullptr, // rotation of union solid - lugShift)); // translation of union solid + plateSolid, // base solid + lugSolid, // union solid + nullptr, // rotation of union solid + lugShift)); // translation of union solid if (verbosityLevel > 3) { G4cout << __PRETTY_FUNCTION__ << " lug " << ithFin << " added with shift (" << lugShift.x() << ", " << lugShift.y() << ", " << lugShift.z() << ")" << G4endl; } @@ -1638,10 +1638,10 @@ namespace mu2e { // add lug fillets if applicable if (plateLugFillet != nullptr) { plateSolid = reg.add(new G4UnionSolid(cacheKey + "_LugFilletUnion_" + std::to_string(ithFin), - plateSolid, // base solid - plateLugFillet, // union solid - stickmanFinRotations.at(ithFin), // rotation of union solid - CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid + plateSolid, // base solid + plateLugFillet, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid } } // end of loop over fins to assemble the plate solid @@ -1657,9 +1657,9 @@ namespace mu2e { int numberOfPlates = tgt->numberOfPlates(); // running z position for plate construction. Start at most negative end and move towards +z double _currentZ = _localCenter.z() - tgt->halfStickmanLength() // end of most negative end - + tgt->stickmanSupportRingLength() // move past support ring at end - + 2 * tgt->spacerHalfLength() ; // move past spacer. - // plate fin, core, and lug is flush on this side + + tgt->stickmanSupportRingLength() // move past support ring at end + + 2 * tgt->spacerHalfLength() ; // move past spacer. + // plate fin, core, and lug is flush on this side if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << " plate Z starts at " << _currentZ << G4endl;} // loop over plates to construct and add them to the target mother. @@ -1674,26 +1674,26 @@ namespace mu2e { if (verbosityLevel > 0) { G4cout << __PRETTY_FUNCTION__ - << " plate " << ithPlate - << " startZ=" << _currentZ - << " centerZ=" << plateCenterZ - << " endZ=" << (_currentZ + tgt->plateLugThickness(ithPlate)) - << " material=" << tgt->plateMaterial(ithPlate) - << G4endl; + << " plate " << ithPlate + << " startZ=" << _currentZ + << " centerZ=" << plateCenterZ + << " endZ=" << (_currentZ + tgt->plateLugThickness(ithPlate)) + << " material=" << tgt->plateMaterial(ithPlate) + << G4endl; } if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ - << " rOut=" << tgt->plateROut(ithPlate) - << " thickness=" << tgt->plateThickness(ithPlate) - << " finWidth=" << tgt->plateFinWidth() - << " finReach=" << tgt->plateFinOuterRadius() - << " lugThickness=" << tgt->plateLugThickness(ithPlate) - << " lug(inner,outer)=(" - << tgt->plateLugInnerRadius() << "," - << tgt->plateLugOuterRadius() << ")" - << " addFilletToPlateCore=" << tgt->addFilletToPlateCore() - << " addFilletToPlateLug=" << tgt->addFilletToPlateLug() - << G4endl; + << " rOut=" << tgt->plateROut(ithPlate) + << " thickness=" << tgt->plateThickness(ithPlate) + << " finWidth=" << tgt->plateFinWidth() + << " finReach=" << tgt->plateFinOuterRadius() + << " lugThickness=" << tgt->plateLugThickness(ithPlate) + << " lug(inner,outer)=(" + << tgt->plateLugInnerRadius() << "," + << tgt->plateLugOuterRadius() << ")" + << " addFilletToPlateCore=" << tgt->addFilletToPlateCore() + << " addFilletToPlateLug=" << tgt->addFilletToPlateLug() + << G4endl; } // get the solid for the plate, using the lambda defined above @@ -1703,14 +1703,14 @@ namespace mu2e { VolumeInfo plateInfo(plateName, currentPlateCenter, prodTargetMotherInfo.centerInWorld); plateInfo.solid = plateSolid; finishNesting(plateInfo // volume info as defined above - ,prodTargetPlateMaterials.at(ithPlate) // material for the plate, can be different for each plate - ,&tgt->productionTargetRotation() // PASSIVE rotation - ,currentPlateCenter // translation to place the plate at the correct z position within the target mother - ,prodTargetMotherInfo.logical // parent volume logical to place the plate in - ,ithPlate // copy number - ,G4Colour::Yellow() // color for visualization - ,"ProductionTarget" // lookup token - ,verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + ,prodTargetPlateMaterials.at(ithPlate) // material for the plate, can be different for each plate + ,&tgt->productionTargetRotation() // PASSIVE rotation + ,currentPlateCenter // translation to place the plate at the correct z position within the target mother + ,prodTargetMotherInfo.logical // parent volume logical to place the plate in + ,ithPlate // copy number + ,G4Colour::Yellow() // color for visualization + ,"ProductionTarget" // lookup token + ,verbosityLevel>1); // pass the verbosity level for printing info in finishNesting // update _currentZ to the END of the next plate, which is the start of the current plate + plateLugThickness(ithPlate) _currentZ += tgt->plateLugThickness(ithPlate); @@ -1718,56 +1718,56 @@ namespace mu2e { // add the rods and spacers G4VSolid* spacerSolid = reg.add(new G4Tubs( - "ProductionTargetSpacerSolid", // - tgt->spacerInnerRadius(), // inner radius - tgt->spacerOuterRadius(), // outer radius - tgt->spacerHalfLength(), // half length - 0., // start angle - CLHEP::twopi)); // sweep angle + "ProductionTargetSpacerSolid", // + tgt->spacerInnerRadius(), // inner radius + tgt->spacerOuterRadius(), // outer radius + tgt->spacerHalfLength(), // half length + 0., // start angle + CLHEP::twopi)); // sweep angle const double rodHalfLength = tgt->rodHalfLength() - stickmanMagicOffset; // make it slightly shorter to avoid precesion issues const double spacerCenterZOffset = tgt->rodHalfLength() - tgt->spacerHalfLength(); for (int ithRod = 0; ithRod < tgt->nStickmanFins(); ++ithRod) { const double rodAngle = tgt->plateFinAngle(ithRod); const CLHEP::Hep3Vector rodCenterInTargetFrame( - tgt->plateCenterToLugCenter() * std::cos(rodAngle), - tgt->plateCenterToLugCenter() * std::sin(rodAngle), - 0.); + tgt->plateCenterToLugCenter() * std::cos(rodAngle), + tgt->plateCenterToLugCenter() * std::sin(rodAngle), + 0.); const CLHEP::Hep3Vector rodCenter = tgt->productionTargetRotation().inverse() * rodCenterInTargetFrame; TubsParams rodParams(0., // rIn - tgt->rodRadius(), // rOut - rodHalfLength); // half length + tgt->rodRadius(), // rOut + rodHalfLength); // half length VolumeInfo rodInfo = nestTubs("ProductionTargetRod_" + std::to_string(ithRod), - rodParams, - prodTargetRodMaterial, - &tgt->productionTargetRotation(), // passive rotation - rodCenter, // translation - prodTargetMotherInfo, // parent volume info for the target mother - ithRod, // copy number - prodTargetVisible, // visibility - G4Colour::Green(), // color for visualization - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); + rodParams, + prodTargetRodMaterial, + &tgt->productionTargetRotation(), // passive rotation + rodCenter, // translation + prodTargetMotherInfo, // parent volume info for the target mother + ithRod, // copy number + prodTargetVisible, // visibility + G4Colour::Green(), // color for visualization + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); if (verbosityLevel > 0) { G4cout << __PRETTY_FUNCTION__ << " rod " << ithRod - << " center=" << rodCenter - << " material=" << tgt->rodMaterial() - << G4endl; + << " center=" << rodCenter + << " material=" << tgt->rodMaterial() + << G4endl; } // place the spacers one spacerHalfLength inward from each rod end const CLHEP::Hep3Vector spacerCenterNegZInTargetFrame( - rodCenterInTargetFrame.x(), - rodCenterInTargetFrame.y(), - -spacerCenterZOffset); + rodCenterInTargetFrame.x(), + rodCenterInTargetFrame.y(), + -spacerCenterZOffset); const CLHEP::Hep3Vector spacerCenterPosZInTargetFrame( - rodCenterInTargetFrame.x(), - rodCenterInTargetFrame.y(), - spacerCenterZOffset); + rodCenterInTargetFrame.x(), + rodCenterInTargetFrame.y(), + spacerCenterZOffset); // transform the spacer centers to the parent frame const CLHEP::Hep3Vector spacerShiftNegZ = tgt->productionTargetRotation().inverse() * spacerCenterNegZInTargetFrame; @@ -1775,13 +1775,13 @@ namespace mu2e { tgt->productionTargetRotation().inverse() * spacerCenterPosZInTargetFrame; if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " spacer " << ithRod - << " target-frame centers neg/pos=" - << spacerCenterNegZInTargetFrame << " / " - << spacerCenterPosZInTargetFrame - << " mother-frame centers neg/pos=" - << spacerShiftNegZ << " / " - << spacerShiftPosZ - << G4endl; + << " target-frame centers neg/pos=" + << spacerCenterNegZInTargetFrame << " / " + << spacerCenterPosZInTargetFrame + << " mother-frame centers neg/pos=" + << spacerShiftNegZ << " / " + << spacerShiftPosZ + << G4endl; } // add the spacers @@ -1789,32 +1789,32 @@ namespace mu2e { VolumeInfo spacerInfoNegZ(spacerNameNegZ, spacerShiftNegZ, prodTargetMotherInfo.centerInWorld); spacerInfoNegZ.solid = spacerSolid; finishNesting(spacerInfoNegZ, // volume info as defined above - prodTargetSpacerMaterial, // material for the spacer - &tgt->productionTargetRotation(), // rotation - spacerShiftNegZ, // translation - prodTargetMotherInfo.logical, // parent logical volume - 2*ithRod, // copy number - G4Colour::Cyan(), // color for visualization - "ProductionTarget", // lookup token - verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + prodTargetSpacerMaterial, // material for the spacer + &tgt->productionTargetRotation(), // rotation + spacerShiftNegZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 2*ithRod, // copy number + G4Colour::Cyan(), // color for visualization + "ProductionTarget", // lookup token + verbosityLevel>1); // pass the verbosity level for printing info in finishNesting std::string spacerNamePosZ = "ProductionTargetSpacerPosZ_" + std::to_string(ithRod); VolumeInfo spacerInfoPosZ(spacerNamePosZ, spacerShiftPosZ, prodTargetMotherInfo.centerInWorld); spacerInfoPosZ.solid = spacerSolid; finishNesting(spacerInfoPosZ, // volume info as defined above - prodTargetSpacerMaterial, // material for the spacer - &tgt->productionTargetRotation(), // rotation - spacerShiftPosZ, // translation - prodTargetMotherInfo.logical, // parent logical volume - 2*ithRod + 1, // copy number - G4Colour::Cyan(), // color for visualization - "ProductionTarget", // lookup token - verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + prodTargetSpacerMaterial, // material for the spacer + &tgt->productionTargetRotation(), // rotation + spacerShiftPosZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 2*ithRod + 1, // copy number + G4Colour::Cyan(), // color for visualization + "ProductionTarget", // lookup token + verbosityLevel>1); // pass the verbosity level for printing info in finishNesting if (verbosityLevel > 0) { G4cout << __PRETTY_FUNCTION__ << " spacer for rod " << ithRod - << " negZ center=" << spacerShiftNegZ - << " posZ center=" << spacerShiftPosZ - << " material=" << tgt->spacerMaterial() - << G4endl; + << " negZ center=" << spacerShiftNegZ + << " posZ center=" << spacerShiftPosZ + << " material=" << tgt->spacerMaterial() + << G4endl; } } // end of loop over rods and spacers @@ -1839,84 +1839,84 @@ namespace mu2e { const double lugStemCenterRadius = tgt->stickmanSupportRingInnerRadius() + lugStemHalfLength; if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " building support ring " << tag - << " positiveEnd=" << positiveEnd - << " ringHalfLength=" << ringHalfLength - << " inner/outer=" << tgt->stickmanSupportRingInnerRadius() - << "/" << tgt->stickmanSupportRingOuterRadius() - << " lugStemHalfLength=" << lugStemHalfLength - << " lugStemCenterRadius=" << lugStemCenterRadius - << G4endl; + << " positiveEnd=" << positiveEnd + << " ringHalfLength=" << ringHalfLength + << " inner/outer=" << tgt->stickmanSupportRingInnerRadius() + << "/" << tgt->stickmanSupportRingOuterRadius() + << " lugStemHalfLength=" << lugStemHalfLength + << " lugStemCenterRadius=" << lugStemCenterRadius + << G4endl; } // ring body without lugs or center hole cutout G4VSolid* supportRingSolid = reg.add(new G4Tubs( - "ProductionTargetSupportRingBase_" + tag, // - 0., // inner radius - tgt->stickmanSupportRingOuterRadius(), // outer radius - ringHalfLength, // half length - 0., // start angle - CLHEP::twopi)); // delta angle + "ProductionTargetSupportRingBase_" + tag, // + 0., // inner radius + tgt->stickmanSupportRingOuterRadius(), // outer radius + ringHalfLength, // half length + 0., // start angle + CLHEP::twopi)); // delta angle // cutout for the center of the ring // such construction avoids potential issue that the lug solid extends within the // inner radius of the ring, which can cause misrepresentation of the geometry G4Tubs* ringCenterCutout = reg.add(new G4Tubs( - "ProductionTargetSupportRingCenterCutout_" + tag, // - 0., // inner radius - tgt->stickmanSupportRingInnerRadius(), // outer radius (avoid overlap with lug stem) - ringHalfLength + stickmanMagicOffset, // half length (add small offset to avoid precision issues) - 0., // start angle - CLHEP::twopi)); // delta angle + "ProductionTargetSupportRingCenterCutout_" + tag, // + 0., // inner radius + tgt->stickmanSupportRingInnerRadius(), // outer radius (avoid overlap with lug stem) + ringHalfLength + stickmanMagicOffset, // half length (add small offset to avoid precision issues) + 0., // start angle + CLHEP::twopi)); // delta angle // lugs at three fin angles. G4Tubs* lugSolid = reg.add(new G4Tubs( - "ProductionTargetSupportRingLug_" + tag, // - 0., // inner radius - tgt->supportRingLugOuterRadius(), // outer radius - ringHalfLength, // half length - 0., // start angle - CLHEP::twopi)); // delta angle + "ProductionTargetSupportRingLug_" + tag, // + 0., // inner radius + tgt->supportRingLugOuterRadius(), // outer radius + ringHalfLength, // half length + 0., // start angle + CLHEP::twopi)); // delta angle // lug stem G4Box* lugStemSolid = reg.add(new G4Box( - "ProductionTargetSupportRingLugStem_" + tag, // - lugStemHalfLength, // half length in x - tgt->supportRingLugOuterRadius(), // half length in y - ringHalfLength)); // half length in z + "ProductionTargetSupportRingLugStem_" + tag, // + lugStemHalfLength, // half length in x + tgt->supportRingLugOuterRadius(), // half length in y + ringHalfLength)); // half length in z for (int ithFin = 0; ithFin < tgt->nStickmanFins(); ++ithFin) { // the lugs are aligned with the fins, so get the fin angle for the current fin to calculate the position of the lugs const double finAngle = tgt->plateFinAngle(ithFin); const CLHEP::Hep3Vector lugStemShift( // rotate the lug stem shift by the fin angle to get the correct position in x and y - lugStemCenterRadius * std::cos(finAngle), - lugStemCenterRadius * std::sin(finAngle), - 0.); + lugStemCenterRadius * std::cos(finAngle), + lugStemCenterRadius * std::sin(finAngle), + 0.); const CLHEP::Hep3Vector lugShift( // rotate the lug shift by the fin angle - tgt->plateCenterToLugCenter() * std::cos(finAngle), - tgt->plateCenterToLugCenter() * std::sin(finAngle), - 0.); + tgt->plateCenterToLugCenter() * std::cos(finAngle), + tgt->plateCenterToLugCenter() * std::sin(finAngle), + 0.); if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " support ring " << tag - << " fin=" << ithFin - << " angle(deg)=" << finAngle/CLHEP::degree - << " lugStemShift=" << lugStemShift - << " lugShift=" << lugShift - << G4endl; + << " fin=" << ithFin + << " angle(deg)=" << finAngle/CLHEP::degree + << " lugStemShift=" << lugStemShift + << " lugShift=" << lugShift + << G4endl; } // first union the lug stem to the ring supportRingSolid = reg.add(new G4UnionSolid( - "ProductionTargetSupportRingLugStemUnion_" + tag + "_" + std::to_string(ithFin), - supportRingSolid, // base solid - lugStemSolid, // union solid - stickmanFinRotations.at(ithFin), // passive rotation - lugStemShift)); // translation + "ProductionTargetSupportRingLugStemUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + lugStemSolid, // union solid + stickmanFinRotations.at(ithFin), // passive rotation + lugStemShift)); // translation // then union the lug to the ring+stem supportRingSolid = reg.add(new G4UnionSolid( - "ProductionTargetSupportRingLugUnion_" + tag + "_" + std::to_string(ithFin), - supportRingSolid, // base solid - lugSolid, // union solid - nullptr, // passive rotation - lugShift)); // translation + "ProductionTargetSupportRingLugUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + lugSolid, // union solid + nullptr, // passive rotation + lugShift)); // translation // add fillet to the lug if applicable // The fillet is added by unioning with an extrusion @@ -1927,14 +1927,14 @@ namespace mu2e { if (tgt->addFilletToSupportRingLug()) { const double filletRadius = tgt->supportRingLugFilletRadius(); const double filletX = std::sqrt( - std::pow(tgt->stickmanSupportRingOuterRadius() + filletRadius, 2) - - std::pow(tgt->supportRingLugOuterRadius() + filletRadius, 2)); + std::pow(tgt->stickmanSupportRingOuterRadius() + filletRadius, 2) + - std::pow(tgt->supportRingLugOuterRadius() + filletRadius, 2)); if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " support ring " << tag - << " fin=" << ithFin - << " filletRadius=" << filletRadius - << " filletX=" << filletX - << G4endl; + << " fin=" << ithFin + << " filletRadius=" << filletRadius + << " filletX=" << filletX + << G4endl; } // get the vertices for the extrusion std::vector filletVertices = { @@ -1945,56 +1945,56 @@ namespace mu2e { // create the extrusion solid for the fillet G4ExtrudedSolid* filletExtrusion = reg.add(new G4ExtrudedSolid( - "ProductionTargetSupportRingLugFilletExtrusion_" + tag + "_" + std::to_string(ithFin), - filletVertices, // vertices - ringHalfLength, // half length in z - G4TwoVector(0., 0.), // shift of polygon center at -half length - 1.0, // scale of polygon at -half length - G4TwoVector(0., 0.), // shift of polygon center at +half length - 1.0)); // scale of polygon at +half length + "ProductionTargetSupportRingLugFilletExtrusion_" + tag + "_" + std::to_string(ithFin), + filletVertices, // vertices + ringHalfLength, // half length in z + G4TwoVector(0., 0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0., 0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length // create the cutout cylinder for the fillet G4Tubs* filletCutout = reg.add(new G4Tubs( - "ProductionTargetSupportRingLugFilletCutout_" + tag + "_" + std::to_string(ithFin), - 0., // inner radius - filletRadius + stickmanMagicOffset, // outer radius (add small offset) - ringHalfLength + stickmanMagicOffset, // half length in z (add small offset) - 0., // start angle - CLHEP::twopi)); // segment angle is full circle + "ProductionTargetSupportRingLugFilletCutout_" + tag + "_" + std::to_string(ithFin), + 0., // inner radius + filletRadius + stickmanMagicOffset, // outer radius (add small offset) + ringHalfLength + stickmanMagicOffset, // half length in z (add small offset) + 0., // start angle + CLHEP::twopi)); // segment angle is full circle // subtract the cutouts from the extrusion G4VSolid* filletUpper = reg.add(new G4SubtractionSolid( - "ProductionTargetSupportRingLugFilletUpper_" + tag + "_" + std::to_string(ithFin), - filletExtrusion, // base solid - filletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(filletX, - tgt->supportRingLugOuterRadius() + filletRadius, - 0.))); // translation of cutout solid + "ProductionTargetSupportRingLugFilletUpper_" + tag + "_" + std::to_string(ithFin), + filletExtrusion, // base solid + filletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(filletX, + tgt->supportRingLugOuterRadius() + filletRadius, + 0.))); // translation of cutout solid G4VSolid* filletSolid = reg.add(new G4SubtractionSolid( - "ProductionTargetSupportRingLugFilletLower_" + tag + "_" + std::to_string(ithFin), - filletUpper, // base solid - filletCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(filletX, - -tgt->supportRingLugOuterRadius() - filletRadius, - 0.))); // translation of cutout solid + "ProductionTargetSupportRingLugFilletLower_" + tag + "_" + std::to_string(ithFin), + filletUpper, // base solid + filletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(filletX, + -tgt->supportRingLugOuterRadius() - filletRadius, + 0.))); // translation of cutout solid // add filletSolid to the union with the support ring supportRingSolid = reg.add(new G4UnionSolid( - "ProductionTargetSupportRingFilletUnion_" + tag + "_" + std::to_string(ithFin), - supportRingSolid, // base solid - filletSolid, // union solid - stickmanFinRotations.at(ithFin), // passive rotation - CLHEP::Hep3Vector(0., 0., 0.))); // no translation + "ProductionTargetSupportRingFilletUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + filletSolid, // union solid + stickmanFinRotations.at(ithFin), // passive rotation + CLHEP::Hep3Vector(0., 0., 0.))); // no translation } // end of fillet construction for the lug } // end of loop over fins to add lugs and fillets to the support ring // add cutout for the center hole supportRingSolid = reg.add(new G4SubtractionSolid( - "ProductionTargetSupportRingLugFilletRingCutout_" + tag, - supportRingSolid, // base solid - ringCenterCutout, // cutout solid - nullptr, // rotation - CLHEP::Hep3Vector(0., 0., 0.))); // translation + "ProductionTargetSupportRingLugFilletRingCutout_" + tag, + supportRingSolid, // base solid + ringCenterCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(0., 0., 0.))); // translation // add cutouts for the through holes if applicable if (tgt->addCutoutToSupportRing()) { @@ -2004,28 +2004,28 @@ namespace mu2e { const double cutoutCenterRadius = 0.5*(tgt->stickmanSupportRingInnerRadius() + tgt->stickmanSupportRingOuterRadius()); // this is z position in the local ring frame const double localCutoutZ = (positiveEnd ? 1. : -1.) * - (tgt->supportRingCutoutOffset() - ringHalfLength + supportRingWallThickness * 0.5 * std::tan(cutoutTiltAngle) ); + (tgt->supportRingCutoutOffset() - ringHalfLength + supportRingWallThickness * 0.5 * std::tan(cutoutTiltAngle) ); // cutout needs to go through the whole wall thickness const double cutoutHalfLength = 0.5 * supportRingWallThickness / std::cos(cutoutTiltAngle) - + tgt->supportRingCutoutInnerRadius() * std::tan(cutoutTiltAngle) - + stickmanMagicOffset; // add small offset to avoid precision issues + + tgt->supportRingCutoutInnerRadius() * std::tan(cutoutTiltAngle) + + stickmanMagicOffset; // add small offset to avoid precision issues if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " support ring " << tag - << " cutoutTilt(deg)=" << cutoutTiltAngle/CLHEP::degree - << " wallThickness=" << supportRingWallThickness - << " cutoutCenterRadius=" << cutoutCenterRadius - << " localCutoutZ=" << localCutoutZ - << " cutoutHalfLength=" << cutoutHalfLength - << G4endl; + << " cutoutTilt(deg)=" << cutoutTiltAngle/CLHEP::degree + << " wallThickness=" << supportRingWallThickness + << " cutoutCenterRadius=" << cutoutCenterRadius + << " localCutoutZ=" << localCutoutZ + << " cutoutHalfLength=" << cutoutHalfLength + << G4endl; } // create the cutout solid for the through hole G4Tubs* supportRingCutout = reg.add(new G4Tubs( - "ProductionTargetSupportRingCutout_" + tag, - 0., // inner radius - tgt->supportRingCutoutInnerRadius(), // outer radius - cutoutHalfLength, // half length in z - 0., // start angle - CLHEP::twopi)); // delta angle + "ProductionTargetSupportRingCutout_" + tag, + 0., // inner radius + tgt->supportRingCutoutInnerRadius(), // outer radius + cutoutHalfLength, // half length in z + 0., // start angle + CLHEP::twopi)); // delta angle const CLHEP::Hep3Vector localZAxis(0., 0., 1.); // loop over the holes and subtract the cutout solid @@ -2033,28 +2033,28 @@ namespace mu2e { const double cutoutPhi = tgt->supportRingCutoutAngle(ithCutout); const CLHEP::Hep3Vector radialAxis(std::cos(cutoutPhi), std::sin(cutoutPhi), 0.); // pointing from ring center radially outward const CLHEP::Hep3Vector cutoutAxis = std::cos(cutoutTiltAngle) * radialAxis // vector pointing outward along the cutout axis - + (positiveEnd ? -1. : 1.) * std::sin(cutoutTiltAngle) * localZAxis; + + (positiveEnd ? -1. : 1.) * std::sin(cutoutTiltAngle) * localZAxis; G4RotationMatrix* cutoutRotation = reg.add(G4RotationMatrix( - CLHEP::HepRotation(cutoutAxis.cross(localZAxis), cutoutAxis.angle(localZAxis)))); // this rotates the cutout from the local z axis to the cutout axis passively + CLHEP::HepRotation(cutoutAxis.cross(localZAxis), cutoutAxis.angle(localZAxis)))); // this rotates the cutout from the local z axis to the cutout axis passively const CLHEP::Hep3Vector cutoutCenter( - cutoutCenterRadius * std::cos(cutoutPhi), - cutoutCenterRadius * std::sin(cutoutPhi), - localCutoutZ); + cutoutCenterRadius * std::cos(cutoutPhi), + cutoutCenterRadius * std::sin(cutoutPhi), + localCutoutZ); if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " support ring " << tag - << " cutout=" << ithCutout - << " phi(deg)=" << cutoutPhi/CLHEP::degree - << " center=" << cutoutCenter - << " axis=" << cutoutAxis - << G4endl; + << " cutout=" << ithCutout + << " phi(deg)=" << cutoutPhi/CLHEP::degree + << " center=" << cutoutCenter + << " axis=" << cutoutAxis + << G4endl; } // make this cutout solid supportRingSolid = reg.add(new G4SubtractionSolid( - "ProductionTargetSupportRingCutoutSubtraction_" + tag + "_" + std::to_string(ithCutout), - supportRingSolid, // base solid - supportRingCutout, // cutout solid - cutoutRotation, // rotation to align the cutout with the cutout axis - cutoutCenter)); // translation + "ProductionTargetSupportRingCutoutSubtraction_" + tag + "_" + std::to_string(ithCutout), + supportRingSolid, // base solid + supportRingCutout, // cutout solid + cutoutRotation, // rotation to align the cutout with the cutout axis + cutoutCenter)); // translation } // end of loop over cutouts for through holes } // end of cutout construction for through holes @@ -2071,56 +2071,56 @@ namespace mu2e { tgt->productionTargetRotation().inverse() * supportRingCenterPosZInTargetFrame; if (verbosityLevel > 2) { G4cout << __PRETTY_FUNCTION__ << " support ring centers target-frame neg/pos=" - << supportRingCenterNegZInTargetFrame << " / " - << supportRingCenterPosZInTargetFrame - << " mother-frame neg/pos=" - << supportRingCenterNegZ << " / " - << supportRingCenterPosZ - << G4endl; + << supportRingCenterNegZInTargetFrame << " / " + << supportRingCenterPosZInTargetFrame + << " mother-frame neg/pos=" + << supportRingCenterNegZ << " / " + << supportRingCenterPosZ + << G4endl; } // add the negative z end ring first VolumeInfo supportRingNegZInfo("ProductionTargetSupportRing_Upstream", - supportRingCenterNegZ, - prodTargetMotherInfo.centerInWorld); + supportRingCenterNegZ, + prodTargetMotherInfo.centerInWorld); supportRingNegZInfo.solid = buildStickmanSupportRingSolid("Upstream", false); finishNesting(supportRingNegZInfo, // volume info as defined above - prodTargetSupportRingMaterial, // material for the support ring - &tgt->productionTargetRotation(), // rotation - supportRingCenterNegZ, // translation - prodTargetMotherInfo.logical, // parent logical volume - 0, // copy number - G4Colour::Green(), - "ProductionTarget", - verbosityLevel>1); + prodTargetSupportRingMaterial, // material for the support ring + &tgt->productionTargetRotation(), // rotation + supportRingCenterNegZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 0, // copy number + G4Colour::Green(), + "ProductionTarget", + verbosityLevel>1); // add the positive z end ring VolumeInfo supportRingPosZInfo("ProductionTargetSupportRing_Downstream", - supportRingCenterPosZ, - prodTargetMotherInfo.centerInWorld); + supportRingCenterPosZ, + prodTargetMotherInfo.centerInWorld); supportRingPosZInfo.solid = buildStickmanSupportRingSolid("Downstream", true); finishNesting(supportRingPosZInfo, // volume info as defined above - prodTargetSupportRingMaterial, // material for the support ring - &tgt->productionTargetRotation(), // rotation - supportRingCenterPosZ, // translation - prodTargetMotherInfo.logical, // parent logical volume - 1, // copy number - G4Colour::Green(), - "ProductionTarget", - verbosityLevel>1); + prodTargetSupportRingMaterial, // material for the support ring + &tgt->productionTargetRotation(), // rotation + supportRingCenterPosZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 1, // copy number + G4Colour::Green(), + "ProductionTarget", + verbosityLevel>1); if (verbosityLevel > 0) { G4cout << __PRETTY_FUNCTION__ << " end rings" - << " negZ center=" << supportRingCenterNegZ - << " posZ center=" << supportRingCenterPosZ - << " material=" << tgt->stickmanSupportRingMaterial() - << G4endl; + << " negZ center=" << supportRingCenterNegZ + << " posZ center=" << supportRingCenterPosZ + << " material=" << tgt->stickmanSupportRingMaterial() + << G4endl; } - //Code following the construction of hayman_v_2_0 to add support structures. The Stickman support wheel - //is the same on design as the Hayman supports, but some errors and inaccuracies in the geometry implementations - //are fixed. + //Code following the construction of hayman_v_2_0 to add support structures. The Stickman support wheel + //is the same on design as the Hayman supports, but some errors and inaccuracies in the geometry implementations + //are fixed. //Keeping it separate in the two branches instead of extracting it to a single helper function may seem redundant - //but is intentional. The Hayman geometry will exist as a benchmark in the future and in principle will not be - //changed or used. So it makes sense to keep it frozen and separated from the Stickman implementation, which may still - //see changes and iterations. Such "duplication" allows the two constructions to exist independently without affecting + //but is intentional. The Hayman geometry will exist as a benchmark in the future and in principle will not be + //changed or used. So it makes sense to keep it frozen and separated from the Stickman implementation, which may still + //see changes and iterations. Such "duplication" allows the two constructions to exist independently without affecting //each other. //Add support structures for the production target @@ -2131,16 +2131,16 @@ namespace mu2e { //create the volume info for the support wheel+rods VolumeInfo suppWheelInfo( "ProductionTargetSupportWheel", localWheelCenter, prodTargetMotherInfo.centerInMu2e()); suppWheelInfo.solid = new G4Tubs("ProductionTargetSupportWheel_wheel", suppWheelParams[0], suppWheelParams[1], - suppWheelParams[2], 0., CLHEP::twopi); - // suppWheelParams, - // suppWheelMaterial, - // 0, - // localWheelCenter, - // prodTargetMotherInfo, - // 0, - // G4Colour::Gray(), - // "PS" - // ); + suppWheelParams[2], 0., CLHEP::twopi); + // suppWheelParams, + // suppWheelMaterial, + // 0, + // localWheelCenter, + // prodTargetMotherInfo, + // 0, + // G4Colour::Gray(), + // "PS" + // ); // add spokes // @@ -2150,9 +2150,9 @@ namespace mu2e { //target info double rTarget = tgt->stickmanSupportRingOuterRadius(); //radius of the support ring to attach to double zTarget = tgt->halfStickmanLength() - - tgt->stickmanSupportRingLength() - + tgt->supportRingCutoutOffset(); - //where along the target to attach, note the asymmetry + - tgt->stickmanSupportRingLength() + + tgt->supportRingCutoutOffset(); + //where along the target to attach, note the asymmetry double smallGap = 0.001; //for adding small offsets to avoid overlaps due to precision //initialize parameter vectors //features on wheel @@ -2186,11 +2186,11 @@ namespace mu2e { : spokeTargetAnglesU[ispoke]*CLHEP::degree; double rWheel = supportWheelRodRadialOffset[ispoke]; // radius of the wheel to attach to CLHEP::Hep3Vector rodCenter(rWheel*std::cos(wheelAngle)/std::cos(targetAngle), - rWheel*std::sin(wheelAngle), - 0.); // rod plane projected to the wheel plane, non-orthogonal cut of a cylinder is an ellipse, - // so need to divide by cos(targetAngle) to get the correct position on the wheel plane + rWheel*std::sin(wheelAngle), + 0.); // rod plane projected to the wheel plane, non-orthogonal cut of a cylinder is an ellipse, + // so need to divide by cos(targetAngle) to get the correct position on the wheel plane const double rodOffset = supportWheelRodOffset[ispoke] - + supportWheelRodPinOffset[ispoke]/std::cos(targetAngle); + + supportWheelRodPinOffset[ispoke]/std::cos(targetAngle); rodCenter += CLHEP::Hep3Vector(std::sin(targetAngle)*rodOffset, 0., std::cos(targetAngle)*rodOffset); if(istream == 0) { //only do once //add the features near the support rods in the bicycle wheel @@ -2203,14 +2203,14 @@ namespace mu2e { CLHEP::Hep3Vector featureCenter(localWheelCenter); //center is wheel center double featureParams[] = {featureRIn, featureROut, tgt->supportWheelHL(), featureAngle - featureArc/2. /*phi0*/, featureArc /*dphi*/}; G4Tubs* featureTubs = new G4Tubs("ProductionTargetSupportFeature_" +std::to_string(ispoke), - featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); + featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelFeature_union_"+std::to_string(ispoke), - suppWheelInfo.solid, featureTubs, 0, featureCenter); + suppWheelInfo.solid, featureTubs, 0, featureCenter); //add the support rod to the wheel G4Tubs* rodTubs = new G4Tubs("ProductionTargetSupportRod_" +std::to_string(ispoke), - 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); + 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelRod_union_"+std::to_string(ispoke), - suppWheelInfo.solid, rodTubs, rodRot, rodCenter); + suppWheelInfo.solid, rodTubs, rodRot, rodCenter); } const int side = (1-2*istream); //+1 or -1 //info about wire connection @@ -2220,8 +2220,8 @@ namespace mu2e { wheelPos += side*supportWheelRodHL[ispoke]*rodAxis; //translate from rod center to edge, this is vector from center to rod center rotated by target angle CLHEP::Hep3Vector rodCenterToWire(std::cos(wheelAngle)*std::cos(targetAngle), - std::sin(wheelAngle), // rorated by target angle around y so y should be unchanged - -std::cos(wheelAngle)*std::sin(targetAngle)); + std::sin(wheelAngle), // rorated by target angle around y so y should be unchanged + -std::cos(wheelAngle)*std::sin(targetAngle)); wheelPos -= supportWheelRodRadius[ispoke]*rodCenterToWire; double zWireRodOffset = (istream == 0) ? supportWheelRodWireOffsetD[ispoke] : supportWheelRodWireOffsetU[ispoke]; wheelPos -= side*zWireRodOffset*rodAxis; @@ -2231,19 +2231,19 @@ namespace mu2e { targetPos = tgt->productionTargetRotation().inverse()*targetPos; //rotate from target frame to mother frame if (verbosityLevel > 1) G4cout << __PRETTY_FUNCTION__ - << "Sanity check for spoke " << ispoke << " on stream " << istream << ":\n" - << "wheel pos " << wheelPos << "\ntarget pos " << targetPos << "\n" - << "spoke length " << (wheelPos-targetPos).mag() << G4endl; + << "Sanity check for spoke " << ispoke << " on stream " << istream << ":\n" + << "wheel pos " << wheelPos << "\ntarget pos " << targetPos << "\n" + << "spoke length " << (wheelPos-targetPos).mag() << G4endl; CLHEP::Hep3Vector spokeAxis((wheelPos-targetPos).unit()); CLHEP::Hep3Vector targetAxis(0.,0.,side); targetAxis = tgt->productionTargetRotation().inverse()*targetAxis; CLHEP::Hep3Vector zax(0.,0.,1.); if(verbosityLevel > 0) G4cout << __PRETTY_FUNCTION__ - << "istream " << istream << " ispoke " << ispoke << G4endl - << "target pos " << targetPos << "\nwheel pos " << wheelPos << G4endl - << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << G4endl - << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << G4endl; + << "istream " << istream << " ispoke " << ispoke << G4endl + << "target pos " << targetPos << "\nwheel pos " << wheelPos << G4endl + << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << G4endl + << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << G4endl; //to remove overlaps where the wire connects, need angle of wire and surface connecting to //remove overlap at target double wireTargetAngle = targetAxis.angle(-1.*spokeAxis); @@ -2251,8 +2251,8 @@ namespace mu2e { targetPos += (deltaLength+0.1)*spokeAxis; //subtract off the length if(verbosityLevel > 0) G4cout << __PRETTY_FUNCTION__ - << "wire target angle " << wireTargetAngle << " delta L " << deltaLength - << " target pos " << targetPos << G4endl; + << "wire target angle " << wireTargetAngle << " delta L " << deltaLength + << " target pos " << targetPos << G4endl; //next remove overlap at rod // @@ -2263,8 +2263,8 @@ namespace mu2e { wheelPos -= (deltaLength+0.1)*spokeAxis; if(verbosityLevel > 0) G4cout << __PRETTY_FUNCTION__ - << "wire rod angle " << wireRodAngle << " delta L " << deltaLength - << " wheel pos " << wheelPos << G4endl; + << "wire rod angle " << wireRodAngle << " delta L " << deltaLength + << " wheel pos " << wheelPos << G4endl; CLHEP::Hep3Vector spokeCenter((wheelPos+targetPos)/2.); double spokeLength = std::abs((wheelPos-targetPos).mag()); @@ -2279,27 +2279,27 @@ namespace mu2e { spokeName << ispoke; VolumeInfo spokeInfo = nestTubs( spokeName.str(), - spokeParams, - spokeMaterial, - spokeRot, - spokeCenter, - prodTargetMotherInfo, - 0, - G4Colour::Gray(), - "PS" - ); + spokeParams, + spokeMaterial, + spokeRot, + spokeCenter, + prodTargetMotherInfo, + 0, + G4Colour::Gray(), + "PS" + ); } //end spokes loop } //end stream loop finishNesting(suppWheelInfo, - suppWheelMaterial, - 0, - localWheelCenter, - prodTargetMotherInfo.logical, - 0, - G4Colour::Gray(), - "PS" - ); + suppWheelMaterial, + 0, + localWheelCenter, + prodTargetMotherInfo.logical, + 0, + G4Colour::Gray(), + "PS" + ); } //end adding support structures } //end constructStickmanTarget diff --git a/Mu2eKinKal/CMakeLists.txt b/Mu2eKinKal/CMakeLists.txt index 51fbad3f3e..3e1ae785bd 100644 --- a/Mu2eKinKal/CMakeLists.txt +++ b/Mu2eKinKal/CMakeLists.txt @@ -8,9 +8,7 @@ cet_make_library( src/KKConstantBField.cc src/KKFitSettings.cc src/KKFitUtilities.cc - src/KKMaterial.cc src/KKSHFlag.cc - src/KKStrawMaterial.cc src/StrawHitUpdaters.cc src/StrawXingUpdater.cc src/WHSIterator.cc diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index c827389d0f..03fc9879a2 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -1,22 +1,9 @@ #include "Offline/TrkReco/fcl/Particle.fcl" BEGIN_PROLOG Mu2eKinKal : { + # # general configuration # - MAT: { - elements : "Offline/TrackerConditions/data/ElementsList.data" - isotopes : "Offline/TrackerConditions/data/IsotopesList.data" - materials : "Offline/TrackerConditions/data/MaterialsList.data" - strawGasMaterialName : "straw-gas" - strawWallMaterialName : "straw-wall" - strawWireMaterialName : "straw-wire" - IPAMaterialName : "HDPE" - STMaterialName : "Target" - IonizationEnergyLossMode : 1 # Moyal mean - SolidScatteringFraction : 0.999999 - GasScatteringFraction : 0.9999999 - ElectronBrehmsFraction : 0.04 - } KKFIT: { PrintLevel : 0 TPOCAPrecision : 1e-4 # mm @@ -446,7 +433,6 @@ Mu2eKinKal : { LHSeedFit : { module_type : LoopHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SampleSurfaces : ["TT_Mid"] @@ -468,7 +454,6 @@ Mu2eKinKal : { LHDriftFit : { module_type : LoopHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings : { @table::Mu2eKinKal.KKFIT SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] # these are additional surfaces; surfaces used in extrapolation are also sampled @@ -493,7 +478,6 @@ Mu2eKinKal : { KLSeedFit : { module_type : KinematicLineFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: @local::Mu2eKinKal.KKFIT FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : { @@ -508,7 +492,6 @@ Mu2eKinKal : { KLDriftFit : { module_type : KinematicLineFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT # DNB I don't know where this time offset difference WRT helical fits comes from. If it's physical, we will need a better @@ -531,7 +514,6 @@ Mu2eKinKal : { CHSeedFit : { module_type : CentralHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : T0 @@ -547,7 +529,6 @@ Mu2eKinKal : { CHDriftFit : { module_type : CentralHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : Detector diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index ab5a4205c7..f46290275d 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -13,7 +13,7 @@ #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/GeometryService.hh" @@ -30,7 +30,7 @@ namespace mu2e { class KKExtrap { public: - explicit KKExtrap(KKExtrapConfig const& exconfig,KKMaterial const& kkmat); + explicit KKExtrap(KKExtrapConfig const& exconfig); // extrapolation functions; these are templated on the type of trajectory template void extrapolate(KKTrack& ktrk) const; template void toTrackerEnds(KKTrack& ktrk) const; @@ -43,18 +43,16 @@ namespace mu2e { private: int debug_; double btol_, intertol_, maxdt_; - KKMaterial const& kkmat_; bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO }; - KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig,KKMaterial const& kkmat) : + KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig) : debug_(extrapconfig.Debug()), btol_(extrapconfig.btol()), intertol_(extrapconfig.interTol()), maxdt_(extrapconfig.MaxDt()), - kkmat_(kkmat), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), @@ -63,8 +61,7 @@ namespace mu2e { template void KKExtrap::extrapolate(KKTrack& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); // define the time direction according to the fit direction inside the tracker auto const& ftraj = ktrk.fitTraj(); @@ -103,10 +100,9 @@ namespace mu2e { template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); - ExtrapolateToZ trackerBack(maxdt_,btol_,kkg.tracker()->back().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerBack(maxdt_,btol_,kkg_h->tracker()->back().center().Z(),debug_); // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't // reflect inside the tracker volume auto const& ftraj = ktrk.fitTraj(); @@ -121,18 +117,18 @@ namespace mu2e { static const SurfaceId tt_back("TT_Back"); // start with the middle - auto midinter = KinKal::intersect(ftraj,kkg.tracker()->middle(),ftraj.range(),intertol_); + auto midinter = KinKal::intersect(ftraj,kkg_h->tracker()->middle(),ftraj.range(),intertol_); if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); if(tofront){ // check the front piece first; that is usually correct // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto frontinter = KinKal::intersect(fhel,kkg.tracker()->front(),fhel.range(),intertol_,fronttdir); + auto frontinter = KinKal::intersect(fhel,kkg_h->tracker()->front(),fhel.range(),intertol_,fronttdir); if(!frontinter.good()){ // start from the middle TimeRange frange = ftraj.range(); if(midinter.good())frange = fronttdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - frontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),frange,intertol_,fronttdir); + frontinter = KinKal::intersect(ftraj,kkg_h->tracker()->front(),frange,intertol_,fronttdir); } if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); } @@ -140,19 +136,19 @@ namespace mu2e { // start from the middle TimeRange brange = ftraj.range(); if(midinter.good())brange = backtdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - auto backinter = KinKal::intersect(ftraj,kkg.tracker()->back(),brange,intertol_,backtdir); + auto backinter = KinKal::intersect(ftraj,kkg_h->tracker()->back(),brange,intertol_,backtdir); if(backinter.good())ktrk.addIntersection(tt_back,backinter); } } template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; // extraplate the fit through the IPA. This will add material effects for each intersection. It will continue till the // track exits the IPA - ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg.DS()->innerProtonAbsorberPtr(),debug_); + ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg_h->DS()->innerProtonAbsorberPtr(),debug_); if(extrapIPA.debug() > 0)std::cout << "extrapolating to IPA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId IPASID("IPA"); @@ -163,8 +159,8 @@ namespace mu2e { if(extrapIPA.intersection().good()){ // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - auto const& IPA = kkg.DS()->innerProtonAbsorberPtr(); - KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_.IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); + auto const& IPA = kkg_h->DS()->innerProtonAbsorberPtr(); + KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_h->IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); if(extrapIPA.debug() > 0){ double dmom, paramomvar, perpmomvar; ipaxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -191,11 +187,11 @@ namespace mu2e { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; // extraplate the fit through the ST. This will add material effects for each foil intersection. It will continue till the // track exits the ST in Z - ExtrapolateST extrapST(maxdt_,btol_,intertol_,*kkg.ST(),debug_); + ExtrapolateST extrapST(maxdt_,btol_,intertol_,*kkg_h->ST(),debug_); auto const& ftraj = ktrk.fitTraj(); double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); auto startdir = ftraj.direction(starttime); @@ -205,7 +201,7 @@ namespace mu2e { if(extrapST.intersection().good()){ // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - KKSTXINGPTR stxingptr = std::make_shared(extrapST.foil(),extrapST.foilId(),*kkmat_.STMaterial(),extrapST.intersection(),reftrajptr,stthick_,extrapST.interTolerance()); + KKSTXINGPTR stxingptr = std::make_shared(extrapST.foil(),extrapST.foilId(),*kkmat_h->STMaterial(),extrapST.intersection(),reftrajptr,stthick_,extrapST.interTolerance()); if(extrapST.debug() > 0){ double dmom, paramomvar, perpmomvar; stxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -231,15 +227,14 @@ namespace mu2e { template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); if(trackerFront.debug() > 0)std::cout << "extrapolating to Tracker " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TrackerSID("TT_Front"); ktrk.extrapolate(tdir,trackerFront); // the last piece appended should cover the necessary range auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto trkfrontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),ktraj.range(),intertol_,tdir); + auto trkfrontinter = KinKal::intersect(ftraj,kkg_h->tracker()->front(),ktraj.range(),intertol_,tdir); if(trkfrontinter.onsurface_){ // dont worry about bounds here ktrk.addIntersection(TrackerSID,trkfrontinter); return true; @@ -249,8 +244,7 @@ namespace mu2e { template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ TSDA(maxdt_,btol_,kkg.DS()->upstreamAbsorber().center().Z(),debug_); + ExtrapolateToZ TSDA(maxdt_,btol_,kkg_h->DS()->upstreamAbsorber().center().Z(),debug_); if(TSDA.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TSDASID("TSDA"); @@ -261,7 +255,7 @@ namespace mu2e { bool retval = epos.Z() < TSDA.zVal(); if(retval){ auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto tsdainter = KinKal::intersect(ftraj,kkg.DS()->upstreamAbsorber(),ktraj.range(),intertol_,tdir); + auto tsdainter = KinKal::intersect(ftraj,kkg_h->DS()->upstreamAbsorber(),ktraj.range(),intertol_,tdir); if(tsdainter.onsurface_)ktrk.addIntersection(TSDASID,tsdainter); } return retval; @@ -269,11 +263,10 @@ namespace mu2e { template void KKExtrap::toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId OPASID("OPA"); TimeRange trange = tdir == TimeDir::forwards ? TimeRange(tstart,ftraj.range().end()) : TimeRange(ftraj.range().begin(),tstart); - auto opainter = KinKal::intersect(ftraj,kkg.DS()->outerProtonAbsorber(),trange,intertol_,tdir); + auto opainter = KinKal::intersect(ftraj,kkg_h->DS()->outerProtonAbsorber(),trange,intertol_,tdir); if(opainter.good()){ ktrk.addIntersection(OPASID,opainter); } diff --git a/Mu2eKinKal/inc/KKFileFinder.hh b/Mu2eKinKal/inc/KKFileFinder.hh deleted file mode 100644 index 67ab289602..0000000000 --- a/Mu2eKinKal/inc/KKFileFinder.hh +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef Mu2eKinKal_KKFileFinder_hh -#define Mu2eKinKal_KKFileFinder_hh - -#include "KinKal/MatEnv/FileFinderInterface.hh" -#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" -#include - -namespace mu2e { - - class KKFileFinder : public MatEnv::FileFinderInterface { - - public: - KKFileFinder(std::string elementsBaseName, - std::string isotopesBaseName, - std::string materialsBaseName): - elementsBaseName_(elementsBaseName), - isotopesBaseName_(isotopesBaseName), - materialsBaseName_(materialsBaseName){} - - std::string matElmDictionaryFileName() const override { - return findFile(elementsBaseName_); - } - - std::string matIsoDictionaryFileName() const override { - return findFile(isotopesBaseName_); - } - - std::string matMtrDictionaryFileName() const override { - return findFile(materialsBaseName_); - } - - std::string findFile( std::string const& path ) const override { - return policy_( path ); - } - - private: - - mutable ConfigFileLookupPolicy policy_; - - std::string elementsBaseName_; - std::string isotopesBaseName_; - std::string materialsBaseName_; - - }; - -} - -#endif /* btrkHelper_KKFileFinder_hh */ diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 3589d7956e..32e47d6aba 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -7,7 +7,7 @@ #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" diff --git a/Mu2eKinKal/inc/KKStrawXing.hh b/Mu2eKinKal/inc/KKStrawXing.hh index 4c09250ea3..1cc2136616 100644 --- a/Mu2eKinKal/inc/KKStrawXing.hh +++ b/Mu2eKinKal/inc/KKStrawXing.hh @@ -6,11 +6,12 @@ #include "KinKal/Detector/ElementXing.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" +#include "Offline/TrackerGeom/inc/Straw.hh" +#include "Offline/DataProducts/inc/StrawId.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" -#include "Offline/DataProducts/inc/StrawId.hh" #include "cetlib_except/exception.h" namespace mu2e { using KinKal::SVEC3; @@ -164,7 +165,7 @@ namespace mu2e { varscale_ = 1.0; } // update the material xings from gas, straw wall, and wire - smat_.findXings(ca_.tpData(),sxconfig_,mxings_); + smat_.findXings(ca_.tpData(),sxconfig_.nsig_,mxings_,sxconfig_.diag_); // update the effect these have on the parameters fparams_ = this->parameterChange(varscale_); } diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index d45316dd45..f6cef066ed 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -102,7 +102,6 @@ namespace mu2e { using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; using Name = fhicl::Name; using Comment = fhicl::Comment; using KinKal::DVEC; @@ -128,7 +127,6 @@ namespace mu2e { fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; // helix module specific config }; @@ -155,7 +153,6 @@ namespace mu2e { int print_; PDGCode::type fpart_; KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix double mass_; // particle mass int PDGcharge_; // PDG particle charge @@ -184,7 +181,6 @@ namespace mu2e { print_(settings().modSettings().printLevel()), fpart_(static_cast(settings().modSettings().fitParticle())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), fixedfield_(false), @@ -246,6 +242,7 @@ namespace mu2e { void CentralHelixFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -319,7 +316,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); // optionally (and if present) add the CaloCluster as a constraint // verify the cluster looks physically reasonable before adding it TODO! Or, let the KKCaloHit updater do it TODO KKCALOHITCOL calohits; @@ -335,7 +332,7 @@ namespace mu2e { // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { // std::cout << "EXTENDING TRACK " << event.id() << " " << index << std::endl; - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); goodfit = goodFit(*kktrk); } diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index c973b93ba9..9f309c0ba8 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -47,7 +47,7 @@ // Mu2eKinKal #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" @@ -103,7 +103,6 @@ namespace mu2e { using KKConfig = Mu2eKinKal::KinKalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; class KinematicLineFit : public art::EDProducer { using Name = fhicl::Name; @@ -133,7 +132,6 @@ namespace mu2e { fhicl::Table mu2eSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; }; @@ -163,7 +161,6 @@ namespace mu2e { PDGCode::type fpart_; TrkFitDirection fdir_; KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix std::array paramconstraints_; double mass_; // particle mass @@ -191,7 +188,6 @@ namespace mu2e { seedmom_(settings().modSettings().seedmom()), fpart_(static_cast(settings().modSettings().fitParticle())), kkfit_(settings().mu2eSettings()), - kkmat_(settings().matSettings()), intertol_(settings().modSettings().interTol()), sampletbuff_(settings().modSettings().sampleTBuff()), sampleinrange_(settings().modSettings().sampleInRange()), @@ -257,13 +253,14 @@ namespace mu2e { kkbf_ = std::make_unique(*bfmgr,*det); // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService GeomHandle kkg_h; - auto const& kkg = *kkg_h; - kkg.surfaces(ssids_,surfacess_to_sample_); + kkg_h->surfaces(ssids_,surfacess_to_sample_); } void KinematicLineFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; + // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -303,7 +300,7 @@ namespace mu2e { KKSTRAWXINGCOL strawxings; strawhits.reserve(strawHitIdxs.size()); strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); //here KKCALOHITCOL calohits; @@ -324,7 +321,7 @@ namespace mu2e { auto kktrk = make_unique(config_,*kkbf_,seedtraj,fpart_,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramconstraints_); auto goodfit = goodFit(*kktrk); if(goodfit && exconfig_.schedule().size() > 0){ - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); } goodfit = goodFit(*kktrk); // extrapolate as required @@ -411,9 +408,9 @@ namespace mu2e { void KinematicLineFit::extrapolate(KKTRK& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; // extrapolate to the extracted CRV. This function should be migrated to KKExtrap TODO - auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg.TCRV(),extrapdebug_); + auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),extrapdebug_); auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TCRVSID("TCRV"); @@ -436,7 +433,7 @@ namespace mu2e { // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); // TODO add DS and shielding material - KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_.STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); + KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); ktrk.addTCRVXing(crvxingptr,tdir); } } while(hadintersection); diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 5e62035125..fca260082a 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -108,7 +108,6 @@ namespace mu2e { using EXINGPTR = std::shared_ptr; using EXINGCOL = std::vector; - using KKMaterialConfig = KKMaterial::Config; using Name = fhicl::Name; using Comment = fhicl::Comment; @@ -125,7 +124,6 @@ namespace mu2e { fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::OptionalTable finalSettings { Name("FinalSettings") }; fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; // LoopHelix module specific config @@ -167,7 +165,6 @@ namespace mu2e { TrkFitDirection fdir_; bool usePDGCharge_; // use the pdg particle charge: otherwise use the helicity and direction to determine the charge KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix double mass_; // particle mass int PDGcharge_; // PDG particle charge @@ -198,7 +195,6 @@ namespace mu2e { useHelixSlope_(settings().slopeSigThreshold(slopeSigThreshold_)), usePDGCharge_(settings().pdgCharge()), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), fixedfield_(false) @@ -224,7 +220,7 @@ namespace mu2e { kkbf_ = std::move(std::make_unique(VEC3(0.0,0.0,bz))); } // setup extrapolation - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); // setup optional fit finalization; this just updates the internals, not the fit result itself if(settings().finalSettings()){ @@ -288,10 +284,9 @@ namespace mu2e { // check the input if(fdir.fitDirection() != TrkFitDirection::FitDirection::downstream && fdir.fitDirection() != TrkFitDirection::FitDirection::upstream) throw cet::exception("RECO") << "mu2e::LoopHelixFit: Unknown helix propagation direction " << fdir.name(); - - // Retrieve event information - // calo geom + // geom GeomHandle calo_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -336,7 +331,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { + if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { if(print_>0) printf("[LoopHelixFit::%s] Failed to create a track\n", __func__); return nullptr; } @@ -358,7 +353,7 @@ namespace mu2e { __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); goodfit = goodFit(*ktrk,seedtraj); // if finalizing, apply that now. if(goodfit && fconfig_.schedule().size() > 0){ diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc index 451c6ecbe7..179331d13c 100644 --- a/Mu2eKinKal/src/RegrowKinematicLine_module.cc +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -54,7 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -79,7 +79,6 @@ namespace mu2e { using Mu2eKinKal::KKFinalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; using SDIS = std::set; using Name = fhicl::Name; @@ -91,7 +90,6 @@ namespace mu2e { fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; }; @@ -130,8 +128,6 @@ namespace mu2e { using DOMAINPTR = std::shared_ptr; using DOMAINCOL = std::set; - using KKMaterialConfig = KKMaterial::Config; - explicit RegrowKinematicLine(const Parameters& settings); void beginRun(art::Run& run) override; void produce(art::Event& event) override; @@ -143,7 +139,6 @@ namespace mu2e { std::unique_ptr kkbf_; Config config_; // refit configuration object, containing the fit schedule KKFIT kkfit_; - KKMaterial kkmat_; art::ProductToken kseedcol_T_; art::ProductToken chcol_T_; art::ProductToken indexmap_T_; @@ -156,7 +151,6 @@ namespace mu2e { debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), kseedcol_T_(consumes(settings().kalSeedCollection())), chcol_T_(consumes(settings().comboHitCollection())), indexmap_T_(consumes(settings().indexMap())), @@ -164,7 +158,7 @@ namespace mu2e { { produces(); produces(); - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); if( fillMCAssns_){ consumes(ksmca_T_); produces (); @@ -185,6 +179,8 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; + // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -218,7 +214,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_.strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 0){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index e4771e5646..a4584dfe52 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -54,7 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -92,7 +92,6 @@ namespace mu2e { fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::Atom extend {Name("Extend"), Comment("Extend the fit") }; @@ -146,7 +145,6 @@ namespace mu2e { std::unique_ptr kkbf_; Config config_; // refit configuration object, containing the fit schedule KKFIT kkfit_; - KKMaterial kkmat_; art::ProductToken kseedcol_T_; art::ProductToken chcol_T_; art::ProductToken cccol_T_; @@ -161,7 +159,6 @@ namespace mu2e { debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), kseedcol_T_(consumes(settings().kalSeedCollection())), chcol_T_(consumes(settings().comboHitCollection())), cccol_T_(mayConsume(settings().caloClusterCollection())), @@ -171,7 +168,7 @@ namespace mu2e { { produces(); produces(); - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); if( fillMCAssns_){ consumes(ksmca_T_); produces (); @@ -192,6 +189,7 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -226,7 +224,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_.strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 1){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; @@ -245,7 +243,7 @@ namespace mu2e { auto ktrk = std::make_unique(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,domains); if(ktrk && ktrk->fitStatus().usable()){ if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl; - if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); + if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); if(ktrk->fitStatus().usable()){ // extrapolate as requested if(extrap_)extrap_->extrapolate(*ktrk); diff --git a/Mu2eKinKal/src/SConscript b/Mu2eKinKal/src/SConscript index 33a145e8e5..1419c781cc 100644 --- a/Mu2eKinKal/src/SConscript +++ b/Mu2eKinKal/src/SConscript @@ -20,6 +20,7 @@ mainlib = helper.make_mainlib([ 'mu2e_GeometryService', 'mu2e_BFieldGeom', 'mu2e_CalorimeterGeom', + 'mu2e_KinKalGeom', 'mu2e_TrackerGeom', 'mu2e_RecoDataProducts', 'mu2e_GlobalConstantsService', @@ -40,7 +41,6 @@ mainlib = helper.make_mainlib([ 'tbb', 'cetlib', 'cetlib_except', - 'mu2e_KinKalGeom', 'KinKal_Fit', 'KinKal_Trajectory', 'KinKal_Geometry', @@ -62,6 +62,7 @@ helper.make_plugins([mainlib, 'mu2e_GeometryService', 'mu2e_BFieldGeom', 'mu2e_CalorimeterGeom', + 'mu2e_KinKalGeom', 'mu2e_TrackerGeom', 'mu2e_RecoDataProducts', 'mu2e_GlobalConstantsService', @@ -87,7 +88,6 @@ helper.make_plugins([mainlib, 'cetlib', 'cetlib_except', 'CLHEP', - 'mu2e_KinKalGeom', 'KinKal_Fit', 'KinKal_Trajectory', 'KinKal_Geometry', diff --git a/ProductionTargetGeom/inc/ProductionTarget.hh b/ProductionTargetGeom/inc/ProductionTarget.hh index c94660563b..0a07a151d0 100644 --- a/ProductionTargetGeom/inc/ProductionTarget.hh +++ b/ProductionTargetGeom/inc/ProductionTarget.hh @@ -78,18 +78,18 @@ namespace mu2e { bool addFilletToPlateCore = false; bool addFilletToPlateLug = false; double plateFilletRadius = 0.0; - + // Support ring fillet parameters bool addFilletToSupportRingLug = false; double supportRingLugFilletRadius = 0.0; - + // Support ring cutout parameters bool addCutoutToSupportRing = false; int nSupportRingCutouts = 0; std::vector supportRingCutoutAngles; double supportRingCutoutInnerRadius = 0.0; double supportRingCutoutTilt = 0.0; - + // Support wheel parameters bool supportsBuild = false; double supportWheelRIn = 0.0; @@ -115,400 +115,400 @@ namespace mu2e { }; class ProductionTarget : virtual public Detector { - public: + public: - // cylinder parameters - std::string tier1TargetType() const {return _tier1TargetType;} - int version() const { return _version; } - double rOut() const { return _rOut; } - double halfLength() const { return _halfLength; } - double envHalfLength() const { return _envelHalfLength; } + // cylinder parameters + std::string tier1TargetType() const {return _tier1TargetType;} + int version() const { return _version; } + double rOut() const { return _rOut; } + double halfLength() const { return _halfLength; } + double envHalfLength() const { return _envelHalfLength; } - // Parameters for version 2 for stiffening parts - double nFins () const { return _nFins; } - double finHeight () const { return _finHeight; } - double finThickness() const { return _finThickness; } - double hubDistUS() const { return _hubDistUS; } - double hubDistDS() const { return _hubDistDS; } - double hubAngleUS() const { return _hubAngleUS; } - double hubAngleDS() const { return _hubAngleDS; } - double hubOverhangUS() const { return _hubOverhangUS;} - double hubOverhangDS() const { return _hubOverhangDS;} - double hubLenUS() const { return _hubLenUS; } - double hubLenDS() const { return _hubLenDS; } - void setHubLenUS(double& aLen ) { _hubLenUS = aLen; } - void setHubLenDS(double& aLen ) { _hubLenDS = aLen; } + // Parameters for version 2 for stiffening parts + double nFins () const { return _nFins; } + double finHeight () const { return _finHeight; } + double finThickness() const { return _finThickness; } + double hubDistUS() const { return _hubDistUS; } + double hubDistDS() const { return _hubDistDS; } + double hubAngleUS() const { return _hubAngleUS; } + double hubAngleDS() const { return _hubAngleDS; } + double hubOverhangUS() const { return _hubOverhangUS;} + double hubOverhangDS() const { return _hubOverhangDS;} + double hubLenUS() const { return _hubLenUS; } + double hubLenDS() const { return _hubLenDS; } + void setHubLenUS(double& aLen ) { _hubLenUS = aLen; } + void setHubLenDS(double& aLen ) { _hubLenDS = aLen; } - // in mu2e coordinates - const CLHEP::Hep3Vector& position() const { return _prodTargetPosition; } + // in mu2e coordinates + const CLHEP::Hep3Vector& position() const { return _prodTargetPosition; } - // this is used to transorm particle momentum and position from - // the PrimaryProtonGun frame to the Mu2e frame - const CLHEP::HepRotation& protonBeamRotation() const { return _protonBeamRotation; } + // this is used to transorm particle momentum and position from + // the PrimaryProtonGun frame to the Mu2e frame + const CLHEP::HepRotation& protonBeamRotation() const { return _protonBeamRotation; } - const CLHEP::Hep3Vector& haymanPosition() const { return _haymanProdTargetPosition; } + const CLHEP::Hep3Vector& haymanPosition() const { return _haymanProdTargetPosition; } - // "passive" rotation, used for placing the production target. - // This is the inverse of protonBeamRotation. - const CLHEP::HepRotation& productionTargetRotation() const { return _protonBeamInverseRotation; } + // "passive" rotation, used for placing the production target. + // This is the inverse of protonBeamRotation. + const CLHEP::HepRotation& productionTargetRotation() const { return _protonBeamInverseRotation; } - const Polycone * getHubsRgtPtr() const {return _pHubsRgtParams/*.get()*/;} - const Polycone * getHubsLftPtr() const {return _pHubsLftParams/*.get()*/;} - const std::map & anchoringPntsRgt() const { return _anchoringPntsRgt; } - const std::map & anchoringPntsLft() const { return _anchoringPntsLft; } + const Polycone * getHubsRgtPtr() const {return _pHubsRgtParams/*.get()*/;} + const Polycone * getHubsLftPtr() const {return _pHubsLftParams/*.get()*/;} + const std::map & anchoringPntsRgt() const { return _anchoringPntsRgt; } + const std::map & anchoringPntsLft() const { return _anchoringPntsLft; } - ~ProductionTarget() override { + ~ProductionTarget() override { - if (!_tier1TargetType.empty()){ - delete _pHubsRgtParams; - delete _pHubsLftParams; + if (!_tier1TargetType.empty()){ + delete _pHubsRgtParams; + delete _pHubsLftParams; + } } - } - ProductionTarget( ProductionTarget const& ) = default; - ProductionTarget( ProductionTarget&& ) = default; - ProductionTarget& operator=(ProductionTarget const& ) = default; - ProductionTarget& operator=(ProductionTarget && ) = default; - - - // - // accessors for hayman v2 - std::string haymanTargetType() const {return _haymanTargetType;} - double halfHaymanLength() const {return _halfHaymanLength;} - std::string targetCoreMaterial() const {return _targetCoreMaterial;} - std::string targetFinMaterial() const {return _targetFinMaterial;} - std::string targetVacuumMaterial() const {return _targetVacuumMaterial;} - std::string supportRingMaterial() const {return _supportRingMaterial;} - std::string spokeMaterial() const {return _spokeMaterial;} - - double rotHaymanX() const {return _rotHaymanX;} - double rotHaymanY() const {return _rotHaymanY;} - double rotHaymanZ() const {return _rotHaymanZ;} - CLHEP::Hep3Vector haymanProdTargetPosition() const {return _haymanProdTargetPosition; } - int numberOfTargetSections() const {return _numberOfTargetSections;} - std::vector startingSectionThickness() const {return _startingSectionThickness;} - std::vector numberOfSegmentsPerSection() const {return _numberOfSegmentsPerSection;} - std::vector thicknessOfSegmentPerSection() const {return _thicknessOfSegmentPerSection;} - std::vector heightOfRectangularGapPerSection() const {return _heightOfRectangularGapPerSection;} - std::vector thicknessOfGapPerSection() const {return _thicknessOfGapPerSection;} - int nHaymanFins() const {return _nHaymanFins;} - std::vector finAngles() const {return _finAngles;} - double haymanFinThickness() const {return _haymanFinThickness;} - double finOuterRadius() const {return _finOuterRadius;} - double supportRingLength() const {return _supportRingLength;} - double supportRingInnerRadius() const {return _supportRingInnerRadius;} - double supportRingOuterRadius() const {return _supportRingOuterRadius;} - double supportRingCutoutThickness() const {return _supportRingCutoutThickness;} - double supportRingCutoutLength() const {return _supportRingCutoutLength;} - - //define the supports for the Hayman_v2 - bool supportsBuild () const {return _supportsBuild ;} - double supportWheelRIn () const {return _supportWheelRIn ;} - double supportWheelROut () const {return _supportWheelROut ;} - double supportWheelHL () const {return _supportWheelHL ;} - std::string supportWheelMaterial() const {return _supportWheelMaterial;} - int nSpokesPerSide () const {return _nSpokesPerSide ;} - //features on the bicycle wheel - const std::vector& supportWheelFeatureAngles() const {return _supportWheelFeatureAngles;} - const std::vector& supportWheelFeatureArcs () const {return _supportWheelFeatureArcs ;} - const std::vector& supportWheelFeatureRIns () const {return _supportWheelFeatureRIns ;} - //rods in the support wheel - const std::vector& supportWheelRodHL () const {return _supportWheelRodHL ;} - const std::vector& supportWheelRodOffset () const {return _supportWheelRodOffset ;} - const std::vector& supportWheelRodPinOffset () const {return _supportWheelRodPinOffset ;} //only used for Stickman_v_1_0 - const std::vector& supportWheelRodRadius () const {return _supportWheelRodRadius ;} - const std::vector& supportWheelRodRadialOffset() const {return _supportWheelRodRadialOffset;} - const std::vector& supportWheelRodWireOffsetD () const {return _supportWheelRodWireOffsetD ;} - const std::vector& supportWheelRodWireOffsetU () const {return _supportWheelRodWireOffsetU ;} - const std::vector& supportWheelRodAngles () const {return _supportWheelRodAngles ;} - - //spoke parameters - const std::vector& spokeTargetAnglesD() const {return _spokeTargetAnglesD;} - const std::vector& spokeTargetAnglesU() const {return _spokeTargetAnglesU;} - double spokeRadius () const {return _spokeRadius ;} - - double productionTargetMotherOuterRadius() const {return _productionTargetMotherOuterRadius;} - double productionTargetMotherHalfLength() const {return _productionTargetMotherHalfLength;} - - //accessors for Stickman_v_1_0 - std::string stickmanTargetType() const {return _stickmanTargetType;} - double halfStickmanLength() const {return _halfStickmanLength;} - int numberOfPlates() const {return _numberOfPlates;} - std::string plateMaterial(int i) const {return _plateMaterial.at(i);} - const std::vector& plateMaterial() const {return _plateMaterial;} - double plateROut(int i) const {return _plateROut.at(i);} - const std::vector& plateROut() const {return _plateROut;} - int nStickmanFins() const {return _nStickmanFins;} - double plateFinAngle(int i) const {return _plateFinAngles.at(i);} - const std::vector& plateFinAngles() const {return _plateFinAngles;} - double plateFinOuterRadius() const {return _plateFinOuterRadius;} - double plateFinWidth() const {return _plateFinWidth;} - double plateCenterToLugCenter() const {return _plateCenterToLugCenter;} - double plateLugInnerRadius() const {return _plateLugInnerRadius;} - double plateLugOuterRadius() const {return _plateLugOuterRadius;} - double plateThickness(int i) const {return _plateThickness.at(i);} - const std::vector& plateThickness() const {return _plateThickness;} - double plateLugThickness(int i) const {return _plateLugThickness.at(i);} - const std::vector& plateLugThickness() const {return _plateLugThickness;} - bool addFilletToPlateCore() const {return _addFilletToPlateCore;} - bool addFilletToPlateLug() const {return _addFilletToPlateLug;} - double plateFilletRadius() const {return _plateFilletRadius;} - std::string rodMaterial() const {return _rodMaterial;} - double rodRadius() const {return _rodRadius;} - double rodHalfLength() const {return _rodHalfLength;} - std::string spacerMaterial() const {return _spacerMaterial;} - double spacerHalfLength() const {return _spacerHalfLength;} - double spacerOuterRadius() const {return _spacerOuterRadius;} - double spacerInnerRadius() const {return _spacerInnerRadius;} - std::string stickmanSupportRingMaterial() const {return _stickmanSupportRingMaterial;} - double stickmanSupportRingLength() const {return _stickmanSupportRingLength;} - double stickmanSupportRingInnerRadius() const {return _stickmanSupportRingInnerRadius;} - double stickmanSupportRingOuterRadius() const {return _stickmanSupportRingOuterRadius;} - double supportRingLugOuterRadius() const {return _supportRingLugOuterRadius;} - bool addFilletToSupportRingLug() const {return _addFilletToSupportRingLug;} - double supportRingLugFilletRadius() const {return _supportRingLugFilletRadius;} - bool addCutoutToSupportRing() const {return _addCutoutToSupportRing;} - int nSupportRingCutouts() const {return _nSupportRingCutouts;} - double supportRingCutoutAngle(int i) const {return _supportRingCutoutAngles.at(i);} - const std::vector& supportRingCutoutAngles() const {return _supportRingCutoutAngles;} - double supportRingCutoutInnerRadius() const {return _supportRingCutoutInnerRadius;} - double supportRingCutoutTilt() const {return _supportRingCutoutTilt;} - double supportRingCutoutOffset() const {return _supportRingCutoutOffset;} - double rotStickmanX() const {return _rotStickmanX;} - double rotStickmanY() const {return _rotStickmanY;} - double rotStickmanZ() const {return _rotStickmanZ;} - CLHEP::Hep3Vector stickmanProdTargetPosition() const {return _stickmanProdTargetPosition;} - - std::string hayman_v_2_0 = "Hayman_v_2_0"; - std::string tier1 = "MDC2018"; - std::string stickman_v_1_0 = "Stickman_v_1_0"; - - CLHEP::Hep3Vector targetPositionByVersion() const { - if (_haymanTargetType == hayman_v_2_0){ - return _haymanProdTargetPosition;} - else if (_tier1TargetType == tier1){ - return _prodTargetPosition;} - else if (_stickmanTargetType == stickman_v_1_0){ - return _stickmanProdTargetPosition;} - else throw cet::exception("BADCONFIG") - << "in ProductionTarget.hh, no valid target specified"<< std::endl; - } - double targetHalfLengthByVersion() const { - if (_haymanTargetType == hayman_v_2_0){ - return _halfHaymanLength;} - else if (_tier1TargetType == tier1){ - return _halfLength;} - else if (_stickmanTargetType == stickman_v_1_0){ - return _halfStickmanLength;} - else throw cet::exception("BADCONFIG") - << "in ProductionTarget.hh, no valid target specified"<< std::endl; - } - - // Configuration method for additional Stickman parameters - void configureStickman(const StickmanConfigParams& configParams); - - //---------------------------------------------------------------- - - private: - friend class ProductionTargetMaker; - - // Private ctr: the class should be only obtained via ProductionTargetFNAL::ProductionTargetMaker. - ProductionTarget(std::string tier1TargetType, int version, double rOut, double halfLength, double rotX, - double rotY, const CLHEP::Hep3Vector& position, - int nFins, - double finHeight, double finThickness, - double hubDistUS, double hubDistDS, - double hubAngleUS, double hubAngleDS, - double hubOverhangUS, double hubOverhangDS ); - - ProductionTarget( - std::string haymanTargetType, int version - ,double productionTargetMotherOuterRadius - ,double productionTargetMotherHalfLength - ,double rOut - ,double halfHaymanLength - ,double rotHaymanX - ,double rotHaymanY - ,double rotHaymanZ - ,const CLHEP::Hep3Vector& haymanProdTargetPosition - ,std::string targetCoreMaterial - ,std::string targetFinMaterial - ,std::string targetVacuumMaterial - ,std::string supportRingMaterial - ,std::string spokeMaterial - ,int numberOfTargetSections - ,std::vector startingSectionThickness - ,std::vector numberOfSegmentsPerSection - ,std::vector thicknessOfSegmentPerSection - ,std::vector heightOfRectangularGapPerSection - ,std::vector thicknessOfGapPerSection - ,int nHaymanFins - ,std::vector finAngles - ,double haymanFinThickness - ,double finOuterRadius - ,double supportRingLength - ,double supportRingInnerRadius - ,double supportRingOuterRadius - ,double supportRingCutoutThickness - ,double supportRingCutoutLength - ); - - ProductionTarget( - std::string stickmanTargetType - ,int version - ,const StickmanEnvelopeParams& envelopeParams - ,const StickmanPlateParams& plateParams - ,const StickmanRodParams& rodParams - ,const StickmanSpacerParams& spacerParams - ,const StickmanSupportRingParams& supportRingParams - ); - - CLHEP::HepRotation _protonBeamRotation; - - // can't return by const ref if invert on the fly so need to store redundant data - CLHEP::HepRotation _protonBeamInverseRotation;// FIXME: should be transient - -// std::unique_ptr _pHubsRgtParams; -// std::unique_ptr _pHubsLftParams; - // Fixme: replace with unique_ptr. - Polycone * _pHubsRgtParams; - Polycone * _pHubsLftParams; - std::map _anchoringPntsRgt; - std::map _anchoringPntsLft; - - CLHEP::Hep3Vector _prodTargetPosition; - - std::string _tier1TargetType; - std::string _haymanTargetType; - - int _version; - double _productionTargetMotherOuterRadius; - double _productionTargetMotherHalfLength; - double _rOut; - double _halfLength; - double _envelHalfLength; - - // version 1+ parameters - int _nFins; - double _finHeight; - double _finThickness; - double _hubDistUS; - double _hubDistDS; - double _hubAngleUS; - double _hubAngleDS; - double _hubOverhangUS; - double _hubOverhangDS; - double _hubLenUS; - double _hubLenDS; - - - // hayman parameters - - - double _halfHaymanLength; - double _rotHaymanX; - double _rotHaymanY; - double _rotHaymanZ; - CLHEP::Hep3Vector _haymanProdTargetPosition; - - std::string _targetCoreMaterial; - std::string _targetFinMaterial; - std::string _targetVacuumMaterial; - std::string _supportRingMaterial; - std::string _spokeMaterial; - int _numberOfTargetSections; - std::vector _startingSectionThickness; - std::vector _numberOfSegmentsPerSection; - std::vector _thicknessOfSegmentPerSection; - std::vector _heightOfRectangularGapPerSection; - std::vector _thicknessOfGapPerSection; - int _nHaymanFins; - std::vector _finAngles; - double _haymanFinThickness; - double _finOuterRadius; - double _supportRingLength; - double _supportRingInnerRadius; - double _supportRingOuterRadius; - double _supportRingCutoutThickness; - double _supportRingCutoutLength; - - //parameters for the support wheel - bool _supportsBuild ; //whether or not to build the supports - double _supportWheelRIn ; - double _supportWheelROut ; - double _supportWheelHL ; //half thickness in z - std::string _supportWheelMaterial; - - //parameters for the non-wheel features near the support rods that are on the wheel - std::vector _supportWheelFeatureAngles; - std::vector _supportWheelFeatureArcs ; - std::vector _supportWheelFeatureRIns ; - - //parameters for rods in the support wheel - std::vector _supportWheelRodHL ; //includes length through the wheel - std::vector _supportWheelRodOffset ; //z offset with respect to the wheel - std::vector _supportWheelRodPinOffset ; //pinhole offset from support wheel center plane, used for Stickman only - std::vector _supportWheelRodRadius ; //radius of the rod - std::vector _supportWheelRodRadialOffset; //radius from the wheel center the rod is centered at - std::vector _supportWheelRodWireOffsetD ; //z offset from the end of the rod the wire connects (downstream) - std::vector _supportWheelRodWireOffsetU ; //z offset from the end of the rod the wire connects (upstream) - std::vector _supportWheelRodAngles ; //angle about the wheel the rod is - - //parameters for the wires (spokes) connecting the support wheel and the target - // std::string _spokeMaterial; (defined above) - int _nSpokesPerSide ; //also constrains the number of rods in the support wheel - std::vector _spokeTargetAnglesD; //angle about the target the wire connects to (downstream) - std::vector _spokeTargetAnglesU; //angle about the target the wire connects to (upstream) - double _spokeRadius ; //radius of the wire - - // Stickman_v_1_0 parameters - std::string _stickmanTargetType; - double _halfStickmanLength; - double _rotStickmanX; - double _rotStickmanY; - double _rotStickmanZ; - CLHEP::Hep3Vector _stickmanProdTargetPosition; - int _numberOfPlates; - std::vector _plateMaterial; - std::vector _plateROut; - int _nStickmanFins; - std::vector _plateFinAngles; - double _plateFinOuterRadius; - double _plateFinWidth; - double _plateCenterToLugCenter; - double _plateLugInnerRadius; - double _plateLugOuterRadius; - std::vector _plateThickness; - std::vector _plateLugThickness; - bool _addFilletToPlateCore; - bool _addFilletToPlateLug; - double _plateFilletRadius; - std::string _rodMaterial; - double _rodRadius; - double _rodHalfLength; - std::string _spacerMaterial; - double _spacerHalfLength; - double _spacerOuterRadius; - double _spacerInnerRadius; - std::string _stickmanSupportRingMaterial; - double _stickmanSupportRingLength; - double _stickmanSupportRingInnerRadius; - double _stickmanSupportRingOuterRadius; - double _supportRingLugOuterRadius; - bool _addFilletToSupportRingLug; - double _supportRingLugFilletRadius; - bool _addCutoutToSupportRing; - int _nSupportRingCutouts; - std::vector _supportRingCutoutAngles; - double _supportRingCutoutInnerRadius; - double _supportRingCutoutTilt; - double _supportRingCutoutOffset; - - // Needed for persistency - template friend class art::Wrapper; - ProductionTarget():_pHubsRgtParams(NULL), _pHubsLftParams(NULL) {} + ProductionTarget( ProductionTarget const& ) = default; + ProductionTarget( ProductionTarget&& ) = default; + ProductionTarget& operator=(ProductionTarget const& ) = default; + ProductionTarget& operator=(ProductionTarget && ) = default; + + + // + // accessors for hayman v2 + std::string haymanTargetType() const {return _haymanTargetType;} + double halfHaymanLength() const {return _halfHaymanLength;} + std::string targetCoreMaterial() const {return _targetCoreMaterial;} + std::string targetFinMaterial() const {return _targetFinMaterial;} + std::string targetVacuumMaterial() const {return _targetVacuumMaterial;} + std::string supportRingMaterial() const {return _supportRingMaterial;} + std::string spokeMaterial() const {return _spokeMaterial;} + + double rotHaymanX() const {return _rotHaymanX;} + double rotHaymanY() const {return _rotHaymanY;} + double rotHaymanZ() const {return _rotHaymanZ;} + CLHEP::Hep3Vector haymanProdTargetPosition() const {return _haymanProdTargetPosition; } + int numberOfTargetSections() const {return _numberOfTargetSections;} + std::vector startingSectionThickness() const {return _startingSectionThickness;} + std::vector numberOfSegmentsPerSection() const {return _numberOfSegmentsPerSection;} + std::vector thicknessOfSegmentPerSection() const {return _thicknessOfSegmentPerSection;} + std::vector heightOfRectangularGapPerSection() const {return _heightOfRectangularGapPerSection;} + std::vector thicknessOfGapPerSection() const {return _thicknessOfGapPerSection;} + int nHaymanFins() const {return _nHaymanFins;} + std::vector finAngles() const {return _finAngles;} + double haymanFinThickness() const {return _haymanFinThickness;} + double finOuterRadius() const {return _finOuterRadius;} + double supportRingLength() const {return _supportRingLength;} + double supportRingInnerRadius() const {return _supportRingInnerRadius;} + double supportRingOuterRadius() const {return _supportRingOuterRadius;} + double supportRingCutoutThickness() const {return _supportRingCutoutThickness;} + double supportRingCutoutLength() const {return _supportRingCutoutLength;} + + //define the supports for the Hayman_v2 + bool supportsBuild () const {return _supportsBuild ;} + double supportWheelRIn () const {return _supportWheelRIn ;} + double supportWheelROut () const {return _supportWheelROut ;} + double supportWheelHL () const {return _supportWheelHL ;} + std::string supportWheelMaterial() const {return _supportWheelMaterial;} + int nSpokesPerSide () const {return _nSpokesPerSide ;} + //features on the bicycle wheel + const std::vector& supportWheelFeatureAngles() const {return _supportWheelFeatureAngles;} + const std::vector& supportWheelFeatureArcs () const {return _supportWheelFeatureArcs ;} + const std::vector& supportWheelFeatureRIns () const {return _supportWheelFeatureRIns ;} + //rods in the support wheel + const std::vector& supportWheelRodHL () const {return _supportWheelRodHL ;} + const std::vector& supportWheelRodOffset () const {return _supportWheelRodOffset ;} + const std::vector& supportWheelRodPinOffset () const {return _supportWheelRodPinOffset ;} //only used for Stickman_v_1_0 + const std::vector& supportWheelRodRadius () const {return _supportWheelRodRadius ;} + const std::vector& supportWheelRodRadialOffset() const {return _supportWheelRodRadialOffset;} + const std::vector& supportWheelRodWireOffsetD () const {return _supportWheelRodWireOffsetD ;} + const std::vector& supportWheelRodWireOffsetU () const {return _supportWheelRodWireOffsetU ;} + const std::vector& supportWheelRodAngles () const {return _supportWheelRodAngles ;} + + //spoke parameters + const std::vector& spokeTargetAnglesD() const {return _spokeTargetAnglesD;} + const std::vector& spokeTargetAnglesU() const {return _spokeTargetAnglesU;} + double spokeRadius () const {return _spokeRadius ;} + + double productionTargetMotherOuterRadius() const {return _productionTargetMotherOuterRadius;} + double productionTargetMotherHalfLength() const {return _productionTargetMotherHalfLength;} + + //accessors for Stickman_v_1_0 + std::string stickmanTargetType() const {return _stickmanTargetType;} + double halfStickmanLength() const {return _halfStickmanLength;} + int numberOfPlates() const {return _numberOfPlates;} + std::string plateMaterial(int i) const {return _plateMaterial.at(i);} + const std::vector& plateMaterial() const {return _plateMaterial;} + double plateROut(int i) const {return _plateROut.at(i);} + const std::vector& plateROut() const {return _plateROut;} + int nStickmanFins() const {return _nStickmanFins;} + double plateFinAngle(int i) const {return _plateFinAngles.at(i);} + const std::vector& plateFinAngles() const {return _plateFinAngles;} + double plateFinOuterRadius() const {return _plateFinOuterRadius;} + double plateFinWidth() const {return _plateFinWidth;} + double plateCenterToLugCenter() const {return _plateCenterToLugCenter;} + double plateLugInnerRadius() const {return _plateLugInnerRadius;} + double plateLugOuterRadius() const {return _plateLugOuterRadius;} + double plateThickness(int i) const {return _plateThickness.at(i);} + const std::vector& plateThickness() const {return _plateThickness;} + double plateLugThickness(int i) const {return _plateLugThickness.at(i);} + const std::vector& plateLugThickness() const {return _plateLugThickness;} + bool addFilletToPlateCore() const {return _addFilletToPlateCore;} + bool addFilletToPlateLug() const {return _addFilletToPlateLug;} + double plateFilletRadius() const {return _plateFilletRadius;} + std::string rodMaterial() const {return _rodMaterial;} + double rodRadius() const {return _rodRadius;} + double rodHalfLength() const {return _rodHalfLength;} + std::string spacerMaterial() const {return _spacerMaterial;} + double spacerHalfLength() const {return _spacerHalfLength;} + double spacerOuterRadius() const {return _spacerOuterRadius;} + double spacerInnerRadius() const {return _spacerInnerRadius;} + std::string stickmanSupportRingMaterial() const {return _stickmanSupportRingMaterial;} + double stickmanSupportRingLength() const {return _stickmanSupportRingLength;} + double stickmanSupportRingInnerRadius() const {return _stickmanSupportRingInnerRadius;} + double stickmanSupportRingOuterRadius() const {return _stickmanSupportRingOuterRadius;} + double supportRingLugOuterRadius() const {return _supportRingLugOuterRadius;} + bool addFilletToSupportRingLug() const {return _addFilletToSupportRingLug;} + double supportRingLugFilletRadius() const {return _supportRingLugFilletRadius;} + bool addCutoutToSupportRing() const {return _addCutoutToSupportRing;} + int nSupportRingCutouts() const {return _nSupportRingCutouts;} + double supportRingCutoutAngle(int i) const {return _supportRingCutoutAngles.at(i);} + const std::vector& supportRingCutoutAngles() const {return _supportRingCutoutAngles;} + double supportRingCutoutInnerRadius() const {return _supportRingCutoutInnerRadius;} + double supportRingCutoutTilt() const {return _supportRingCutoutTilt;} + double supportRingCutoutOffset() const {return _supportRingCutoutOffset;} + double rotStickmanX() const {return _rotStickmanX;} + double rotStickmanY() const {return _rotStickmanY;} + double rotStickmanZ() const {return _rotStickmanZ;} + CLHEP::Hep3Vector stickmanProdTargetPosition() const {return _stickmanProdTargetPosition;} + + std::string hayman_v_2_0 = "Hayman_v_2_0"; + std::string tier1 = "MDC2018"; + std::string stickman_v_1_0 = "Stickman_v_1_0"; + + CLHEP::Hep3Vector targetPositionByVersion() const { + if (_haymanTargetType == hayman_v_2_0){ + return _haymanProdTargetPosition;} + else if (_tier1TargetType == tier1){ + return _prodTargetPosition;} + else if (_stickmanTargetType == stickman_v_1_0){ + return _stickmanProdTargetPosition;} + else throw cet::exception("BADCONFIG") + << "in ProductionTarget.hh, no valid target specified"<< std::endl; + } + double targetHalfLengthByVersion() const { + if (_haymanTargetType == hayman_v_2_0){ + return _halfHaymanLength;} + else if (_tier1TargetType == tier1){ + return _halfLength;} + else if (_stickmanTargetType == stickman_v_1_0){ + return _halfStickmanLength;} + else throw cet::exception("BADCONFIG") + << "in ProductionTarget.hh, no valid target specified"<< std::endl; + } + + // Configuration method for additional Stickman parameters + void configureStickman(const StickmanConfigParams& configParams); + + //---------------------------------------------------------------- + + private: + friend class ProductionTargetMaker; + + // Private ctr: the class should be only obtained via ProductionTargetFNAL::ProductionTargetMaker. + ProductionTarget(std::string tier1TargetType, int version, double rOut, double halfLength, double rotX, + double rotY, const CLHEP::Hep3Vector& position, + int nFins, + double finHeight, double finThickness, + double hubDistUS, double hubDistDS, + double hubAngleUS, double hubAngleDS, + double hubOverhangUS, double hubOverhangDS ); + + ProductionTarget( + std::string haymanTargetType, int version + ,double productionTargetMotherOuterRadius + ,double productionTargetMotherHalfLength + ,double rOut + ,double halfHaymanLength + ,double rotHaymanX + ,double rotHaymanY + ,double rotHaymanZ + ,const CLHEP::Hep3Vector& haymanProdTargetPosition + ,std::string targetCoreMaterial + ,std::string targetFinMaterial + ,std::string targetVacuumMaterial + ,std::string supportRingMaterial + ,std::string spokeMaterial + ,int numberOfTargetSections + ,std::vector startingSectionThickness + ,std::vector numberOfSegmentsPerSection + ,std::vector thicknessOfSegmentPerSection + ,std::vector heightOfRectangularGapPerSection + ,std::vector thicknessOfGapPerSection + ,int nHaymanFins + ,std::vector finAngles + ,double haymanFinThickness + ,double finOuterRadius + ,double supportRingLength + ,double supportRingInnerRadius + ,double supportRingOuterRadius + ,double supportRingCutoutThickness + ,double supportRingCutoutLength + ); + + ProductionTarget( + std::string stickmanTargetType + ,int version + ,const StickmanEnvelopeParams& envelopeParams + ,const StickmanPlateParams& plateParams + ,const StickmanRodParams& rodParams + ,const StickmanSpacerParams& spacerParams + ,const StickmanSupportRingParams& supportRingParams + ); + + CLHEP::HepRotation _protonBeamRotation; + + // can't return by const ref if invert on the fly so need to store redundant data + CLHEP::HepRotation _protonBeamInverseRotation;// FIXME: should be transient + + // std::unique_ptr _pHubsRgtParams; + // std::unique_ptr _pHubsLftParams; + // Fixme: replace with unique_ptr. + Polycone * _pHubsRgtParams; + Polycone * _pHubsLftParams; + std::map _anchoringPntsRgt; + std::map _anchoringPntsLft; + + CLHEP::Hep3Vector _prodTargetPosition; + + std::string _tier1TargetType; + std::string _haymanTargetType; + + int _version; + double _productionTargetMotherOuterRadius; + double _productionTargetMotherHalfLength; + double _rOut; + double _halfLength; + double _envelHalfLength; + + // version 1+ parameters + int _nFins; + double _finHeight; + double _finThickness; + double _hubDistUS; + double _hubDistDS; + double _hubAngleUS; + double _hubAngleDS; + double _hubOverhangUS; + double _hubOverhangDS; + double _hubLenUS; + double _hubLenDS; + + + // hayman parameters + + + double _halfHaymanLength; + double _rotHaymanX; + double _rotHaymanY; + double _rotHaymanZ; + CLHEP::Hep3Vector _haymanProdTargetPosition; + + std::string _targetCoreMaterial; + std::string _targetFinMaterial; + std::string _targetVacuumMaterial; + std::string _supportRingMaterial; + std::string _spokeMaterial; + int _numberOfTargetSections; + std::vector _startingSectionThickness; + std::vector _numberOfSegmentsPerSection; + std::vector _thicknessOfSegmentPerSection; + std::vector _heightOfRectangularGapPerSection; + std::vector _thicknessOfGapPerSection; + int _nHaymanFins; + std::vector _finAngles; + double _haymanFinThickness; + double _finOuterRadius; + double _supportRingLength; + double _supportRingInnerRadius; + double _supportRingOuterRadius; + double _supportRingCutoutThickness; + double _supportRingCutoutLength; + + //parameters for the support wheel + bool _supportsBuild ; //whether or not to build the supports + double _supportWheelRIn ; + double _supportWheelROut ; + double _supportWheelHL ; //half thickness in z + std::string _supportWheelMaterial; + + //parameters for the non-wheel features near the support rods that are on the wheel + std::vector _supportWheelFeatureAngles; + std::vector _supportWheelFeatureArcs ; + std::vector _supportWheelFeatureRIns ; + + //parameters for rods in the support wheel + std::vector _supportWheelRodHL ; //includes length through the wheel + std::vector _supportWheelRodOffset ; //z offset with respect to the wheel + std::vector _supportWheelRodPinOffset ; //pinhole offset from support wheel center plane, used for Stickman only + std::vector _supportWheelRodRadius ; //radius of the rod + std::vector _supportWheelRodRadialOffset; //radius from the wheel center the rod is centered at + std::vector _supportWheelRodWireOffsetD ; //z offset from the end of the rod the wire connects (downstream) + std::vector _supportWheelRodWireOffsetU ; //z offset from the end of the rod the wire connects (upstream) + std::vector _supportWheelRodAngles ; //angle about the wheel the rod is + + //parameters for the wires (spokes) connecting the support wheel and the target + // std::string _spokeMaterial; (defined above) + int _nSpokesPerSide ; //also constrains the number of rods in the support wheel + std::vector _spokeTargetAnglesD; //angle about the target the wire connects to (downstream) + std::vector _spokeTargetAnglesU; //angle about the target the wire connects to (upstream) + double _spokeRadius ; //radius of the wire + + // Stickman_v_1_0 parameters + std::string _stickmanTargetType; + double _halfStickmanLength; + double _rotStickmanX; + double _rotStickmanY; + double _rotStickmanZ; + CLHEP::Hep3Vector _stickmanProdTargetPosition; + int _numberOfPlates; + std::vector _plateMaterial; + std::vector _plateROut; + int _nStickmanFins; + std::vector _plateFinAngles; + double _plateFinOuterRadius; + double _plateFinWidth; + double _plateCenterToLugCenter; + double _plateLugInnerRadius; + double _plateLugOuterRadius; + std::vector _plateThickness; + std::vector _plateLugThickness; + bool _addFilletToPlateCore; + bool _addFilletToPlateLug; + double _plateFilletRadius; + std::string _rodMaterial; + double _rodRadius; + double _rodHalfLength; + std::string _spacerMaterial; + double _spacerHalfLength; + double _spacerOuterRadius; + double _spacerInnerRadius; + std::string _stickmanSupportRingMaterial; + double _stickmanSupportRingLength; + double _stickmanSupportRingInnerRadius; + double _stickmanSupportRingOuterRadius; + double _supportRingLugOuterRadius; + bool _addFilletToSupportRingLug; + double _supportRingLugFilletRadius; + bool _addCutoutToSupportRing; + int _nSupportRingCutouts; + std::vector _supportRingCutoutAngles; + double _supportRingCutoutInnerRadius; + double _supportRingCutoutTilt; + double _supportRingCutoutOffset; + + // Needed for persistency + template friend class art::Wrapper; + ProductionTarget():_pHubsRgtParams(NULL), _pHubsLftParams(NULL) {} }; } diff --git a/RecoDataProducts/inc/TrkStraw.hh b/RecoDataProducts/inc/TrkStraw.hh index 3af028030b..306f2447e2 100644 --- a/RecoDataProducts/inc/TrkStraw.hh +++ b/RecoDataProducts/inc/TrkStraw.hh @@ -8,7 +8,8 @@ #define RecoDataProducts_TrkStraw_HH #include "Offline/DataProducts/inc/StrawId.hh" #include "Offline/RecoDataProducts/inc/StrawFlag.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" +#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" #include "KinKal/Trajectory/ClosestApproachData.hh" namespace mu2e { @@ -25,7 +26,7 @@ namespace mu2e { _dmom(dmom) { double wallpath, gaspath, wirepath; - _pcalc = smat.pathLengths(pocadata,caconfig,wallpath,gaspath,wirepath); + _pcalc = smat.pathLengths(pocadata,caconfig.nsig_,wallpath,gaspath,wirepath); _wallpath = wallpath; _gaspath = gaspath; _wirepath = wirepath; diff --git a/TrackerConditions/data/MaterialsList.data b/TrackerConditions/data/MaterialsList.data index 9d2b64e7e6..01b9b5451b 100644 --- a/TrackerConditions/data/MaterialsList.data +++ b/TrackerConditions/data/MaterialsList.data @@ -1,5 +1,5 @@ # -# This is an exhaustive list of materials for entry into the Conditions DB +# This is an exhaustive list of materials used in Mu2e KinKal # # Do not remove the following line it's a tag ! #Materials_list @@ -29,4 +29,9 @@ cathode 5.2 0. 0. +2 86.4e-2 Copper 0 13.6e-2 straw-wall 1 -10 # effective straw wall + gas material, assuming uniform mass distribution straw_15um 1.927e-2 0.0 0.0 +2 91.3e-2 straw-wall 1 8.7e-2 straw-gas 1 -10 -20 -30 20.0 1.0 solid straw_25um 3.068e-2 0.0 0.0 +2 94.6e-2 straw-wall 1 5.4e-2 straw-gas 1 -10 -20 -30 20.0 1.0 solid +# DS materials +NbTi 6.5 0.0 0.0 +2 0.65 Niobium 0.35 Titanium 0 -10 -20 -30 20.0 1.0 solid +DS1Coil 3.35 0.0 0.0 +2 0.67 Aluminum 0.33 NbTi 0 -10 -20 -30 20.0 1.0 solid +DS2Coil 3.031 0.0 0.0 +2 0.813 Aluminum 0.187 NbTi 0 -10 -20 -30 20.0 1.0 solid +DSSStainless 8.02 0.0 0.0 +5 0.02 Manganese 0.01 Silicon 0.19 Chromium 0.1 Nickel 0.68 Iron 0 -10 -20 -30 20.0 1.0 solid diff --git a/TrackerConditions/fcl/prolog.fcl b/TrackerConditions/fcl/prolog.fcl index edc76089bd..2fb1239208 100644 --- a/TrackerConditions/fcl/prolog.fcl +++ b/TrackerConditions/fcl/prolog.fcl @@ -1,18 +1,4 @@ BEGIN_PROLOG -Mu2eMaterial : { - verbose : 0 - # Location of dictionary files for the material model. - elements : "Offline/TrackerConditions/data/ElementsList.data" - isotopes : "Offline/TrackerConditions/data/IsotopesList.data" - materials : "Offline/TrackerConditions/data/MaterialsList.data" - strawGasMaterialName : "straw-gas" - strawWallMaterialName : "straw-wall" - strawWireMaterialName : "straw-wire" - dahlLynchScatteringFraction : 0.995 - intersectionTolerance : 0.001 - strawElementOffset : 0.25 - maximumIntersectionRadiusFraction : 0.96 -} Mu2eDetector : { verbose : 0 diff --git a/fcl/standardServices.fcl b/fcl/standardServices.fcl index bbf5dcde2f..bea908f6c4 100644 --- a/fcl/standardServices.fcl +++ b/fcl/standardServices.fcl @@ -5,6 +5,7 @@ #include "Offline/fcl/minimalMessageService.fcl" #include "Offline/DbService/fcl/prolog.fcl" #include "Offline/ProditionsService/fcl/prolog.fcl" +#include "Offline/KinKalGeom/fcl/prolog.fcl" BEGIN_PROLOG @@ -37,11 +38,12 @@ Services : { # define services for specific tasks: these are components needed for # a complete job Core : { - message : @local::default_message - GeometryService : { - inputFile : "Offline/Mu2eG4/geom/geom_common.txt" - bFieldFile : "Offline/Mu2eG4/geom/bfgeom_v01.txt" - simulatedDetector : { tool_type: "Mu2e" } + message : @local::default_message + GeometryService : { + inputFile : "Offline/Mu2eG4/geom/geom_common.txt" + bFieldFile : "Offline/Mu2eG4/geom/bfgeom_v01.txt" + simulatedDetector : { tool_type: "Mu2e" } + KinKalMaterial : @local::KinKalGeom.KKMaterial } GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } DbService : @local::DbEmpty