-
Notifications
You must be signed in to change notification settings - Fork 97
Calo Surface Steps #1757
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Draft
sophiemiddleton
wants to merge
15
commits into
Mu2e:main
Choose a base branch
from
sophiemiddleton:interdet
base: main
Could not load branches
Branch not found: {{ refName }}
Loading
Could not load tags
Nothing to show
Loading
Are you sure you want to change the base?
Some commits from the old base branch may be removed from the timeline,
and old review comments may become outdated.
Draft
Calo Surface Steps #1757
Changes from all commits
Commits
Show all changes
15 commits
Select commit
Hold shift + click to select a range
cee3459
initial development
sophieMu2e 9162348
added new Kinkal vols
sophieMu2e 3b9699e
added many surfaces, need positions
sophieMu2e 8219beb
updated coords
sophieMu2e 20c7dc8
updated surfaces
sophieMu2e 0a5a985
builds
sophieMu2e 10d9916
compiles
sophieMu2e 5525394
prints out, not tested
sophieMu2e a5e9665
added surface steps from front and back of calo
sophieMu2e 7c5bb42
added
sophieMu2e 556b4cf
works
sophieMu2e 1d47bee
autoplay
sophieMu2e 23ab8d1
updated Calo
sophieMu2e 959d8f6
Merge branch 'main' of github.com:Mu2e/Offline into interdet
sophieMu2e 2610df8
removed
sophieMu2e File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,48 @@ | ||
| // | ||
| // Define the nominal calorimeter boundary and reference surfaces, used to extrapolate and sample KinKal track fits, and to build | ||
| // the passive materials in the fit | ||
| // original author: Sophie Middleton (2025) | ||
| // | ||
| #ifndef KinKalGeom_Calo_hh | ||
| #define KinKalGeom_Calo_hh | ||
| #include "KinKal/Geometry/Cylinder.hh" | ||
| #include "KinKal/Geometry/Disk.hh" | ||
| #include "KinKal/Geometry/Annulus.hh" | ||
| #include <memory> | ||
| namespace mu2e { | ||
| namespace KinKalGeom { | ||
| class Calo { | ||
| public: | ||
| using CylPtr = std::shared_ptr<KinKal::Cylinder>; | ||
| using DiskPtr = std::shared_ptr<KinKal::Disk>; | ||
| // default constructor with nominal geometry | ||
| Calo(); | ||
| // return by reference | ||
| auto const& EMC_Disk_0_Outer() const { return *EMC_Disk_0_Outer_;} | ||
| auto const& EMC_Disk_0_Inner() const { return *EMC_Disk_0_Inner_;} | ||
| auto const& EMC_Disk_1_Inner() const { return *EMC_Disk_1_Inner_;} | ||
| auto const& EMC_Disk_1_Outer() const { return *EMC_Disk_1_Outer_;} | ||
|
|
||
| auto const& EMC_Disk_0_Front() const { return *EMC_Disk_0_Front_;} | ||
| auto const& EMC_Disk_1_Front() const { return *EMC_Disk_1_Front_;} | ||
| auto const& EMC_Disk_0_Back() const { return *EMC_Disk_0_Back_;} | ||
| auto const& EMC_Disk_1_Back() const { return *EMC_Disk_1_Back_;} | ||
|
|
||
| // return by ptr | ||
| auto const& EMC_Disk_0_OuterPtr() const { return EMC_Disk_0_Outer_;} | ||
| auto const& EMC_Disk_0_InnerPtr() const { return EMC_Disk_0_Inner_;} | ||
| auto const& EMC_Disk_1_InnerPtr() const { return EMC_Disk_1_Inner_;} | ||
| auto const& EMC_Disk_1_OuterPtr() const { return EMC_Disk_1_Outer_;} | ||
| auto const& EMC_Disk_0_FrontPtr() const { return EMC_Disk_0_Front_;} | ||
| auto const& EMC_Disk_1_FrontPtr() const { return EMC_Disk_1_Front_;} | ||
| auto const& EMC_Disk_0_BackPtr() const { return EMC_Disk_0_Back_;} | ||
| auto const& EMC_Disk_1_BackPtr() const { return EMC_Disk_1_Back_;} | ||
|
|
||
| private: | ||
| CylPtr EMC_Disk_0_Inner_, EMC_Disk_0_Outer_ , EMC_Disk_1_Inner_, EMC_Disk_1_Outer_; | ||
| DiskPtr EMC_Disk_0_Front_, EMC_Disk_0_Back_, EMC_Disk_1_Front_, EMC_Disk_1_Back_; | ||
| }; | ||
| } | ||
| } | ||
|
|
||
| #endif |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,31 @@ | ||
| #include "Offline/KinKalGeom/inc/Calo.hh" | ||
| namespace mu2e { | ||
| namespace KinKalGeom { | ||
| using KinKal::VEC3; | ||
| using KinKal::Cylinder; | ||
| using KinKal::Disk; | ||
|
|
||
| // D0 = 11842 -10175.0 = 1667 | ||
| // D1 = 13220 -10175.0 = 3045 | ||
| /* | ||
| double calorimeter.caloDiskRadiusIn = 335; | ||
| double calorimeter.caloDiskRadiusOut = 719; | ||
| half length = 192.295 | ||
|
|
||
|
|
||
| Cylinder(VEC3 const& axis, VEC3 const& center, double radius, double halflen ); | ||
| Disk(VEC3 const& norm,VEC3 const& udir, VEC3 const& center, double radius) | ||
| */ | ||
| Calo::Calo() : | ||
| EMC_Disk_0_Inner_ { std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),335,192.295)}, | ||
| EMC_Disk_0_Outer_ { std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, | ||
| EMC_Disk_1_Inner_ { std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),335,192.295)}, | ||
| EMC_Disk_1_Outer_ { std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, | ||
|
|
||
| EMC_Disk_0_Front_ { std::make_shared<Disk>(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667-192.295),719.)}, | ||
| EMC_Disk_0_Back_ { std::make_shared<Disk>(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+192.295),719.)}, | ||
| EMC_Disk_1_Front_ { std::make_shared<Disk>(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045-192.295),719.)}, | ||
| EMC_Disk_1_Back_ { std::make_shared<Disk>(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045+192.295 ),719.)} | ||
| {} | ||
| } | ||
| } | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,160 @@ | ||
| // predicate to extrapolate to calo | ||
| // Sophie Middleton (2025) | ||
| #ifndef Mu2eKinKal_ExtrapolateCalo_hh | ||
| #define Mu2eKinKal_ExtrapolateCalo_hh | ||
| #include "Offline/Mu2eKinKal/inc/KKTrack.hh" | ||
| #include "KinKal/General/TimeDir.hh" | ||
| #include "KinKal/General/TimeRange.hh" | ||
| #include "KinKal/Geometry/Intersection.hh" | ||
| #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" | ||
| #include "Offline/KinKalGeom/inc/Calo.hh" | ||
| #include "cetlib_except/exception.h" | ||
| #include <memory> | ||
| #include <vector> | ||
| #include <limits> | ||
| namespace mu2e { | ||
| using KinKal::TimeDir; | ||
| using KinKal::TimeRange; | ||
| using KinKalGeom::Calo; | ||
| using KinKal::Annulus; | ||
| using KinKal::Intersection; | ||
| using AnnPtr = std::shared_ptr<KinKal::Annulus>; | ||
| using CaloDisk = std::vector<AnnPtr>; | ||
| using CylPtr = std::shared_ptr<KinKal::Cylinder>; | ||
| class ExtrapolateCalo { | ||
| public: | ||
| ExtrapolateCalo() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10), | ||
| d0zmin_(std::numeric_limits<float>::max()), | ||
| d0zmax_(-std::numeric_limits<float>::max()), | ||
| d1zmin_(std::numeric_limits<float>::max()), | ||
| d1zmax_(-std::numeric_limits<float>::max()), | ||
| rmin_(std::numeric_limits<float>::max()), | ||
| rmax_(-std::numeric_limits<float>::max()), | ||
| debug_(0){} | ||
|
|
||
| ExtrapolateCalo(double maxdt, double dptol, double intertol, Calo const& calo, int debug=0) : | ||
| maxDt_(maxdt), dptol_(dptol), intertol_(intertol), | ||
| d0zmin_( (calo.EMC_Disk_0_Outer().center() - calo.EMC_Disk_0_Outer().axis()*calo.EMC_Disk_0_Outer().halfLength()).Z()), | ||
| d0zmax_( (calo.EMC_Disk_0_Outer().center() + calo.EMC_Disk_0_Outer().axis()*calo.EMC_Disk_0_Outer().halfLength()).Z()), | ||
| d1zmin_( (calo.EMC_Disk_1_Outer().center() - calo.EMC_Disk_1_Outer().axis()*calo.EMC_Disk_1_Outer().halfLength()).Z()), | ||
| d1zmax_( (calo.EMC_Disk_1_Outer().center() + calo.EMC_Disk_1_Outer().axis()*calo.EMC_Disk_1_Outer().halfLength()).Z()), | ||
| rmin_( calo.EMC_Disk_0_Inner().radius()), rmax_(calo.EMC_Disk_0_Outer().radius()), | ||
| disks_(2),debug_(debug) {} | ||
| // interface for extrapolation | ||
| double maxDt() const { return maxDt_; } | ||
| double dpTolerance() const { return dptol_; } | ||
| double interTolerance() const { return intertol_; } | ||
| double d0zmin() const { return d0zmin_; } | ||
| double d0zmax() const { return d0zmax_; } | ||
| double d1zmin() const { return d1zmin_; } | ||
| double d1zmax() const { return d1zmax_; } | ||
| double rmin() const { return rmin_; } | ||
| double rmax() const { return rmax_; } | ||
| auto const& disks () const { return disks_; } | ||
| auto const& intersection() const { return inter_; } | ||
|
|
||
| int debug() const { return debug_; } | ||
| // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time | ||
| template <class KTRAJ> bool needsExtrapolation(KinKal::ParticleTrajectory<KTRAJ> const& fittraj, TimeDir tdir) const; | ||
| // reset between tracks | ||
| void reset() const { inter_ = Intersection(); sid_ = SurfaceId(); ann_ = AnnPtr();} | ||
| // find the nearest disk to a z positionin a given z direction | ||
| size_t nearestDisk(double zpos, double zdir) const; | ||
| private: | ||
| double maxDt_; // maximum extrapolation time | ||
| double dptol_; // fractional momentum tolerance | ||
| double intertol_; // intersection tolerance (mm) | ||
| mutable Intersection inter_; // cache of most recent intersection | ||
| mutable SurfaceId sid_; // cache of most recent intersection disk SID | ||
| mutable AnnPtr ann_; // cache of most recent intersection disk surface | ||
| // cache of front and back Z positions | ||
| double d0zmin_, d0zmax_; // z range of disk0 volume | ||
| double d1zmin_, d1zmax_; // z range of disk1 volume | ||
| double rmin_, rmax_; // inner and outer radii of the anuli | ||
| CaloDisk disks_; // disks | ||
| int debug_; // debug level | ||
| }; | ||
|
|
||
| template <class KTRAJ> bool ExtrapolateCalo::needsExtrapolation(KinKal::ParticleTrajectory<KTRAJ> const& fittraj, TimeDir tdir) const { | ||
| // we are answering the question: did the segment last added to this extrapolated track hit a calo disk or not? If so we are done | ||
| // extrapolating (for now) and we want to find all the intersections in that piece. If not, and if we're still inside or heading towards the | ||
| // disks, keep going. | ||
| auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); | ||
| // add a small buffer to the test range to prevent re-intersection with the same piece | ||
| static const double epsilon(1e-7); // small difference to avoid re-intersecting | ||
| if(ktraj.range().range() <= epsilon) return true; // keep going if the step is very small | ||
| auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; | ||
| auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); | ||
| auto vel = ktraj.velocity(stime); // physical velocity | ||
| if(tdir == TimeDir::backwards) vel *= -1.0; | ||
| auto spos = ktraj.position3(stime); | ||
| auto epos = ktraj.position3(etime); | ||
| if(debug_ > 2)std::cout << "calo extrap tdir " << tdir << " start z " << spos.Z() << " end z " << epos.Z() << " zvel " << vel.Z() << " rho " << spos.Rho() << std::endl; | ||
| // stop if the particle is heading away from the calo | ||
| if( (vel.Z() > 0 && spos.Z() > d1zmax_ ) || (vel.Z() < 0 && spos.Z() < d0zmin_)){ | ||
| reset(); // clear any cache | ||
| if(debug_ > 1)std::cout << "Heading away from calo: done" << std::endl; | ||
| return false; | ||
| } | ||
| // if the particle is going in the right direction but haven't yet reached the calo in Z just keep going | ||
| if( (vel.Z() > 0 && epos.Z() < d0zmin_) || (vel.Z() < 0 && epos.Z() > d1zmax_) ){ | ||
| reset(); | ||
| if(debug_ > 2)std::cout << "Heading towards calo, z " << spos.Z()<< std::endl; | ||
| return true; | ||
| } | ||
| // if we get to here we are in the correct Z range. Test disks. | ||
| int idisk = nearestDisk(spos.Z(),vel.Z()); | ||
| if(idisk >= (int)disks_.size())return true; | ||
| if(debug_ > 2)std::cout << "Looping on disks " << std::endl; | ||
| int ddisk = vel.Z() > 0.0 ? 1 : -1; // iteration direction | ||
| auto trange = tdir == TimeDir::forwards ? TimeRange(stime,ktraj.range().end()) : TimeRange(ktraj.range().begin(),stime); | ||
| // loop over disks in the z range of this piece | ||
| while(idisk >= 0 && idisk < (int)disks_.size() && (disks_[idisk]->center().Z() - epos.Z())*ddisk < 0.0){ | ||
| auto diskptr = disks_[idisk]; | ||
| if(debug_ > 2)std::cout << "disk " << idisk << " z " << diskptr->center().Z() << std::endl; | ||
| auto newinter = KinKal::intersect(ktraj,*diskptr,trange,intertol_,tdir); | ||
| if(debug_ > 2)std::cout << "calo disk inter " << newinter << std::endl; | ||
| if(newinter.good()){ | ||
| // update the cache | ||
| inter_ = newinter; | ||
| ann_ = disks_[idisk]; | ||
| //sid_ = SurfaceId(SurfaceIdEnum::calo_Foils,idisk); //FIXME | ||
| if(debug_ > 0)std::cout << "Good calo disk " << newinter << " sid " << sid_ << std::endl; | ||
| return false; | ||
| } | ||
| idisk += ddisk; // otherwise continue loopin on disks | ||
| } | ||
| // no more intersections: keep extending in Z till we clear the calo | ||
| reset(); | ||
| if(debug_ > 1)std::cout << "Extrapolating to calo edge, z " << spos.Z() << std::endl; | ||
| if(vel.Z() > 0.0) | ||
| return spos.Z() < d1zmax_; | ||
| else | ||
| return spos.Z() > d0zmin_; | ||
| } | ||
|
|
||
| size_t ExtrapolateCalo::nearestDisk(double zpos, double zvel) const { | ||
| size_t retval = disks_.size(); | ||
| if(zvel > 0.0){ // going forwards in z | ||
| for(auto idisk= disks_.begin(); idisk != disks_.end(); idisk++){ | ||
| auto const& diskptr = *idisk; | ||
| if(diskptr->center().Z() > zpos){ | ||
| retval = std::distance(disks_.begin(),idisk); | ||
| break; | ||
| } | ||
| } | ||
| } else { | ||
| for(auto idisk= disks_.rbegin(); idisk != disks_.rend(); idisk++){ | ||
| auto const& diskptr = *idisk; | ||
| if(diskptr->center().Z() < zpos){ | ||
| auto jdisk = idisk.base()-1; // points to the equivalent forwards object | ||
| retval = std::distance(disks_.begin(),jdisk); | ||
| break; | ||
| } | ||
| } | ||
| } | ||
| return retval; | ||
| } | ||
|
|
||
| } | ||
| #endif |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
z position here is what seems to be incorrect - need to understand why