search

Mu2e / Offline

Offline software for the Mu2e experiment
Apache License 2.0
8 stars 82 forks source link

Kalman filter adaptation of cosmic fit #66

Closed ryuwd closed 4 years ago

ryuwd commented 4 years ago

Kalman filter adaptation of @sophiemiddleton 's cosmic fit

Opening this issue from the project since it looks like this is now being developed / discussed.

Needed for #61 .

ryuwd commented 4 years ago

@sophiemiddleton anything I can do to help with this? I have pulled your latest changes from CosmicDriftUpdate.

sophiemiddleton commented 4 years ago

Well I think the main thing we need to do is have a set way of paramerizing the tracks using the TrkParam and Traj style, like its done for HelixTraj. I had an exchange wit Dave earlier and we have some ideas as to how to proceed. I still dont know if the TrkStrawHit is needed. I have some code to do this (its commented out in the latest push of CosmicTrackfit) it works but seg faults if I comment back in a line in the vprop in BTrk/TrkStrawHit constructor. This is what has prevent ed the use of this infrastructure. If you found a way around that it would be useful but tbh may be impossible so if you have a spare bit of time feel free to have a think but dont put it before more important work

brownd1978 commented 4 years ago

Any problem with vprop can be debugged. TrkStrawHit has interfaces to all the physics corrections plus Kalman fit eventually straw alignment, it's working and maintained in other (helix) context, I would like to reuse it here.

Kalman fit for cosmic tracks is a big project. If multiple people are going to work on it we need to have a common design. I have requirements also coming from compatibility with the 6-parameter kinematic helix fit. I will be away the first 2 weeks of Dec., we could have a brief discussion this week then something more concrete in 3 weeks. Dave

On Mon, Nov 25, 2019 at 11:13 AM Sophie Middleton notifications@github.com wrote:

Well I think the main thing we need to do is have a set way of paramerizing the tracks using the TrkParam and Traj style, like its done for HelixTraj. I had an exchange wit Dave earlier and we have some ideas as to how to proceed. I still dont know if the TrkStrawHit is needed. I have some code to do this (its commented out in the latest push of CosmicTrackfit) it works but seg faults if I comment back in a line in the vprop in BTrk/TrkStrawHit constructor. This is what has prevent ed the use of this infrastructure. If you found a way around that it would be useful but tbh may be impossible so if you have a spare bit of time feel free to have a think but dont put it before more important work

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/Mu2e/Offline/issues/66?email_source=notifications&email_token=ABAH57622E4AIPL7SILOLZLQVQPUZA5CNFSM4JRMPXOKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFDPOPY#issuecomment-558298943, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABAH572EXYSOAFSGX3657ZDQVQPUZANCNFSM4JRMPXOA .

-- David Nathan Brown Dave_Brown@lbl.gov Office Phone (510) 486-7261 Fax 495-2957 Lawrence Berkeley National Lab MS 50R5008 (50-6026C) Berkeley, CA 94720

sophiemiddleton commented 4 years ago

Hi Dave,

Yes a discussion might be good- we can draft out what we need to do. I think I should do most of the Cosmic side of things (obviously with your advice) and Ryun can help me once we have a plan, if he wants. I would prefer to keep it just a few people working on it so the 2 (Me and Dave) of us should be enough but if guided assistance can be given, Ryun could help.

The issue with vprop is probably something to do with the StrawResponse not being passed correctly or something, if I comment out that line the rest of the code seems to work fine, I think with your help we could easily debug that and I have the code to use the TrkStrawHit once we get that issue fixed so the upgrade from that point shouldnt be too hard.

I am open to a quick online chat about this sometime before you go away. I am fairly free so whenever suits you would be good for me. Let me know, Sophie

ryuwd commented 4 years ago

Well I think the main thing we need to do is have a set way of paramerizing the tracks using the TrkParam and Traj style, like its done for HelixTraj. I had an exchange wit Dave earlier and we have some ideas as to how to proceed. I still dont know if the TrkStrawHit is needed. I have some code to do this (its commented out in the latest push of CosmicTrackfit) it works but seg faults if I comment back in a line in the vprop in BTrk/TrkStrawHit constructor. This is what has prevent ed the use of this infrastructure. If you found a way around that it would be useful but tbh may be impossible so if you have a spare bit of time feel free to have a think but dont put it before more important work

I can try running it through GDB. Is it this line: https://github.com/sophiemiddleton/Offline/commit/9783e56d7973eb6cfabaf8ac5d80ecf2b1e607f8#r36128263 ?

ryuwd commented 4 years ago

I think the problem is here: https://github.com/sophiemiddleton/Offline/commit/9783e56d7973eb6cfabaf8ac5d80ecf2b1e607f8#r36128758

ryuwd commented 4 years ago

I've made a class that allows configuration of Millepede-II and writing of the binary file using the Mille class provided by Millepede. The binary file is the input to the main solver application called ./pede. I've drafted also a Plane alignment module which 'registers' all the planes as alignable rigid bodies with 6 alignment parameters (translation vector + 3 rotation angles). Also very easy to extend to panels.

Now just need to feed it the residuals (from kalman inverse filter step) and residual derivatives with respect to the alignment parameters and track parameters, so I can test. ./pede will also require another text file for applying constraints etc... which can be made further down the line.

The code is on my GitHub fork here.

I tried to also get the MLE re-fitting stuff working, but it's just too slow and there's no way I'm going to get it to converge 100% of the time.

sophiemiddleton commented 4 years ago

excellent good work- I plan to get going on the Kalman infrastructure changes over the weekend (currently bank holiday today and yesterday so been a slow week). I'll let you know how things are changing by early next week. I think the overall changes will take some time though

ryuwd commented 4 years ago

excellent good work- I plan to get going on the Kalman infrastructure changes over the weekend (currently bank holiday today and yesterday so been a slow week). I'll let you know how things are changing by early next week. I think the overall changes will take some time though

OK, great! Did you get around to writing up what you and Dave discussed? I'd be very happy to assist with helping write the code / parts of the infrastructure if this is possible, and if the design decisions have been made. I realise it's a big-ish project, and given the direct dependency of the Cosmic Alignment on this I'd like to help as much as possible.

Sorry - I realise now it's your bank holiday, but are there any small-ish tasks I could start with now / work on over the weekend? (Lectures are off due to UCU strikes (until 6th) and I've finished almost all of my coursework workload for this semester, so I have a lot more time than usual to dedicate to Mu2e.) Cheers

brownd1978 commented 4 years ago

Hi Ryun and Sophie, Sophie and I agreed on a provisional parameterization of the KF linear trajectory class. The (geometric) components will be: -theta and phi of the direction vector WRT to the tracker coordinate system -signed distance of closest approach (DOCA) to the tracker system origin

The trajectory parameterization, the tracker geometry, and the residual definition provides enough information to write the code to compute the derivatives. Each residual is the signed DOCA between the local trajectory, minus the predicted drift distance. The drift distance is a function of the drift time, and in this pass we will ignore time-domain effects (which are 2nd order), so you can consider it fixed. Thus the residual derivatives are just the change in DOCA WRT the alignment parameters: dDOCA/dPi, where Pi are (say) the 6 degrees of freedom for the plane alignment. These can be computed as a function of the track parameters (above) and the position and direction of the wire. I suggest first writing a function to compute DOCA (you can use the TwoLinePCA class in Mu2eUtilitites or equivalent) given this parameterization, and a utility to compute the derivatives numerically by changing the alignment parameters. Those derivatives can be used directly in the alignment, or better, used to check an eventual analytic calculation of the derivatives. You can test the calculations with a standalone program based on the Mu2e geometry.

I'll be on vacation with minimal connectivity through Dec. 13. I'll try to keep up with email. Hopefully this is enough to keep everybody busy. cheers, Dave

On Fri, Nov 29, 2019 at 8:06 AM ryuwd notifications@github.com wrote:

excellent good work- I plan to get going on the Kalman infrastructure changes over the weekend (currently bank holiday today and yesterday so been a slow week). I'll let you know how things are changing by early next week. I think the overall changes will take some time though

OK, great! Did you get around to writing up what you and Dave discussed? I'd be very happy to assist with helping write the code / parts of the infrastructure if this is possible, and if the design decisions have been made. I realise it's a big-ish project, and given the direct dependency of the Cosmic Alignment on this I'd like to help as much as possible.

Sorry - I realise now it's your bank holiday, but are there any small-ish tasks I could start with now / work on over the weekend? (Lectures are off due to UCU strikes (until 6th) and I've finished almost all of my coursework workload for this semester, so I have a lot more time than usual to dedicate to Mu2e.)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/Mu2e/Offline/issues/66?email_source=notifications&email_token=ABAH575VHWACSJETBL46YQDQWE4XDA5CNFSM4JRMPXOKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFPGCKA#issuecomment-559833384, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABAH577GYMUFUIZM3T5Y3ZTQWE4XDANCNFSM4JRMPXOA .

-- David Nathan Brown Dave_Brown@lbl.gov Office Phone (510) 486-7261 Fax 495-2957 Lawrence Berkeley National Lab MS 50R5008 (50-6026C) Berkeley, CA 94720

sophiemiddleton commented 4 years ago

Yeah sounds good. I think Dave mostly summarised what we discussed, but this is what I have in my raw notes:

T0 → interesting points, currently we fixed t0 in the minuit fit. If ignore t0 we can go about in terms of a geometric fit, if not we need to upgrade to kinematic fit

The Cosmic BTrk interface:

-We need position and direction -Direction  theta and phi

Momentum: -TrkMomCalcuator used by KalRep -Get momentum but needs some B field, even with field off, just set this to very small value -momentum is just external constant in field off scenario

Curvature:

-Can set to return a constant value -Field Mu2eG4 common/current  DS field -option number 2 -set so that momentum is O(GeV) -EL and MCS  neg. for cosmics

DifAlgebra:

-HelixTraj model -HelixParams come from DifParams

KalFit Base this new fitter on top of patrec (current) fitter…

Not that detailed but I think Dave's last paragraph helps you with what you are wanting to do over the next week or so.

Meanwhile, as I begin implementing things if a nice task pops out that you could assist with I'll let you know.

Cheers, Sophie

ryuwd commented 4 years ago

Hi Dave, Sophie,

Thank you both very much for the detailed replies. I'd been thinking about the analytical solution for the derivatives the last few days and I've had a go at calculating the residual/DOCA partial derivatives analytically with sympy in this python script. The script generates C code which can be compiled in Offline. The alignment parameters and track, wire info can then be used to calculate the residual derivatives analytically (can't test until I crosscheck with numerical results).

I will adapt it to work with the straight track traj parameters as well, also will get started on the numerical calculations using the geometry (which I will cross check with the maths / functions my script have generated). Cheers

sophiemiddleton commented 4 years ago

Looks good, I'm going to need a week I think to get the kalman stuff up and working, I have some dummy code but its just lengthy process joining all the loose ends.....will let you know when its at a reasonable standard

sophiemiddleton commented 4 years ago

This seems to be a lot more code than I realized to be able to replicate what the Helix Traj/Params does. I have a few meetings today and for the rest of the week but I hope to at least have all the place holders ready by the end of the week. It will probably be most of next week then checking it all compiles and making sure its doing what we think its doing. I am aiming to have something like KalFit_module for the CosmicTracks. This will consume CosmicTrackSeedCollection and produce CosmicKalCollection. I think the situation here is simpler than the Helix fit so maybe we dont need all the infrastructure that is needed there. I have currently added: CosmicTraj, CosmicParams, (in my Offline but similar to HelixTraj/Params), CosmicVal, CosmicCov (offline verisons), CosmicTrkUtilities, (used to make some track stuff) this is all for the trajectory fitting. I dont know about KalSegment (in Offline ) I have added in an accessor to a cosmic and comsic_cov but Im not sure that will be used. We also need RecoProducts like CosmicKalSeed, CosmicKalFitData. I need to add some overloads to KalRep, but this will be my own version of that, there's no need to re-do the whole code, I think just init and another function. I will have to think how to access the TrkMomCalculator. I know we have a plan for what to do there yet but I havent put it into practice yet (thats probably todays job, in between meetings).

This is what I've done so far. There are math things I need to calculate for example and I need to understand some of the propagation in the Traj. I need also to link all these things together, i'm sure there will be some redundancy here. I need to figure out what we need to call the "centre" of the track and things like that (mostly reinterpretation). The actual CosmicKalFit helper function will be based on KalFit, with something like MakeTrack, but there's a lot of other functionality in there which I need to decide if we need.

Anyway, just updating on where things are at present.

Cheers, Sophie

ryuwd commented 4 years ago

Great! Thanks for this.

I've been digging around on the subject of max likelihood estimation with biased models i.e. not modelling detector misalignment in our case. I think we can take the Variance-Covariance matrix of the errors (for each straw hit) and calculate the Cholesky decomposition L.

I've read a few resources that state that L can then be used to calculate uncorrelated, unbiased estimates of the residuals by rotating the biased residuals for a track using the inverse of L.

I think this might be a good simple way to first order to calculate unbiased residuals with the current Minuit drift fit implementation (without re-fitting), in the meantime whilst the Kalman fit is developed. Though I definitely think of course in the long term it will still be better and necessary to use a Kalman fit to find cosmic track trajectories.

I'm now able to run Millepede on output from my module, and am trying to gather enough track information to start being able to calculate alignment parameters for planes. (I can't get it to actually accept my data yet..!) I need to also work on figuring out how many tracks Millepede will need to be able to solve for alignment parameters. I think it could be possibly quite a lot. I have about 4000 so far from my last run (just on my computer though). I maybe should start running my module on the grid soon.

sophiemiddleton commented 4 years ago

Well, all sounds like good progress. Can you forward me the resources so I can have a read too. Thanks.

Yeah I think the grid is a good option- do you know how to run on the virtual machines? All the info is on the Mu2e wiki for grid running but sometimes there's things you miss the first time round so let me know if you have any issues.

sophiemiddleton commented 4 years ago

The transition to Kalman might be a little longer. I need to think some more about the track parameterization. All the infrastrucutre is ready (ish) but its more just ensuring that we get the best representation of the track.