Feature #454
Review beam maker code
100%
Description
Just want to remind myself to have a look
Related issues
Updated by Rogers, Chris about 12 years ago
Relevant files:
MapPyBeamMaker
test_MapPyBeamMaker
MapPyCreateSpill
test_MapPyCreateSpill
Updated by Rogers, Chris about 12 years ago
- Assignee changed from Rogers, Chris to Carlisle, Timothy
General¶
A good first iteration. I like the fact that you split the process to create the number of events and to create the actual distribution into two mappers. Of course, I have some comments and bugs.
Comments¶
Please have a look at the required comment syntax at
http://micewww.pp.rl.ac.uk/projects/maus/wiki/Coding_style
your comments are all far too short. You should define for each parameter, what the parameter does and what types it takes, any return value, plus any in-function logic.
ErrorHandler¶
- You should handle errors using the methods in ErrorHandler class
http://micewww.pp.rl.ac.uk/embedded/maus/doxygen_framework/html/classcore_1_1ErrorHandler_1_1ErrorHandler.html
This allows us to set at run time whether errors go into json document, go to user, or something else. If usage is unclear, let me know and I will explain/improve the comments.
Testing¶
- I don't quite feel that you've got the difference between tests that we run to make sure the code is correct (which happend in test_MapPyCreateSpill, etc) and tests where we check that the user input is correct (which we do in the actual runtime code, i.e. MapPyCreateSpill)
MapPyCreateSpill¶
Comments¶
- Class description in the comment like """---""" is not good enough! Please comment properly here.
- Presumably birth() does nothing as the comment is just """Birth"""? Even if it does nothing you should say so.
- etc, comments clearly need to be improved
Code¶
MapPyCreateSpill.process()¶
- In process you have
print " "
statement. Shouldn't be there. - when you make an error, the correct format should be spill['errors']['MapPyCreateSpill']['error message']
- BUG indentation error:
if 'errors' not in spill: spill['errors'] = {} spill['errors']['pre-existing_branch'] = 'pre_existing MCbranch' return json.dumps(spill)
should be
if 'errors' not in spill: spill['errors'] = {} spill['errors']['pre-existing_branch'] = 'pre_existing MCbranch' return json.dumps(spill)
this stuff would all go away if you used the error handler class.
test_MapPyCreateSpill¶
- please remove the print statement from test_Nparticles; but I don't really understand what you are trying to test with this function - surely better test would be to test that
self.n_particles == config_doc["CreateSpill_Nparticles"]
If you want to check CreateSpill_Nparticles >= 0, this would be better to go in the birth() function?
- you really really should check that the number of particles created is equal to "CreateSpill_Nparticles"; you should probably see what happens when n particles is 0, negative, 1, >1 and check code behaves as expected - these seem to be the obvious boundary cases to me.
MapPyBeamMaker¶
Comments
- as above. Please fix.
Code¶
MapPyBeamMaker.process()¶
- BUG
b_z = 4.# Tesla kappa = 0.15 * ( (b_z) / ( self._central_pz ) ) # m^-1 beta = 1. / math.fabs( kappa ) gamma = 2. / beta # print "beta = " + str(beta*100) + " mm" mass = 105.658367 # MeV/c ## defined centrally?
particle["particle_id"] = 13# PDG PID
This stuff should never be hard-coded. Please make them into datacards (except mass, see below).
Nb: If you import xboa.Common, there is a system of units (xboa.Common.units), some physical constants (xboa.Common.constants) and pre-defined masses, pdg pid system, string naming system for many particle types (xboa.Common.pdg_pid_to_mass, xboa.Common.pdg_pid_to_name). Note the bug http://micewww.pp.rl.ac.uk/issues/453.- BUG
_xmax = 5. * math.sqrt( s_xx ) _pmax = 5. * math.sqrt( s_pp )
This doesn't give the correct emittance, right? - BUG unit momentum is not a unit vector
MapPyBeamMaker.death()¶
- BUG Should return true
test_MapPyBeamMaker¶
- BUG you never tested that MapPyBeamMaker.process() actually does anything! The main point of these tests is to check that MapPyBeamMaker.process() creates particles with x,y,z,px,py,pz,t,energy, and that they are:
(a) physical
(b) with correct distribution
but you never checked this.
- test_BeamMakerCards() - what is the point of this test? We don't plan to run tests on the datacards like this - again, this sort of test should be (and is) in the MapPyBeamMaker.birth() process
Updated by Rogers, Chris about 12 years ago
Nb: Tim, could you have a look at this stuff and then bounce the issue back to me (make Assigned to Chris Rogers) when you are finished. I'll try to ring and chat later anyway.
Updated by Carlisle, Timothy about 12 years ago
Sure thing.
I rang earlier, but will try again after 2.30pm.
Updated by Carlisle, Timothy about 12 years ago
Take Two has now been pushed...hopefully the tests are better!
A few questions Chris:
1) ErrorHandler:
I'm unsure how to check if exceptions have been raised by the ErrorHandler
see: MapPyBeamMaker.py L135 - L142, where the exceptions are sent to the handler
& test_MapPyBeamMaker.py L54-57 & L228-230, where I wish to check if an exception was raised (the old style error flags are there at the moment).
also an error is flagged in the ErrorHandler code during test_MapPyCreateSpill.py - I'm probably doing something wrong here!
2) The 5sigma range L164/165 in BeamMaker.py was lifted from the G4MICE equivalent - won't the hit & miss method catch any 'bad' particles that are outside of the phase space?
Further comments/corrections welcome,
Thanks
Updated by Tunnell, Christopher about 12 years ago
I agree with Tim here: this is a hit/miss MC. It pulls from a uniform distribution over +/- 5 sigma then sees if that value lies under the multivariate gaussian curve. Maybe it requires some documentation.
There are other ways to generate random gaussian numbers that are gauranteed to terminate and more computationally efficient. This will give the correct answer I believe for a multivariate gaussian beam though.
Updated by Rogers, Chris almost 12 years ago
- Assignee changed from Carlisle, Timothy to Rogers, Chris
Sorry - I should have looked at this before, didn't see it because was still assigned to Tim
Updated by Rogers, Chris almost 12 years ago
What is the name of your branch? So I can check it out
Updated by Tunnell, Christopher almost 12 years ago
said it for him so he knew how to complete the request ;)
Updated by Tunnell, Christopher almost 12 years ago
Oh, I found out something interesting today. So currently the logic for the beammaker says that it is doing a hit miss MC (ie. sample uniformly and if it's not under the curve, toss it out). Does it mention that in the documentation? I forget... anyways.
There is actually a better way to do this than a hit miss MC. You can transform a multivariate gaussian into a space where you can just use a random 1D number using 'box muller transformations' (still figuring that out).
http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
anyways, so I wrote something that let's numpy do a lot of the work and I think it's a lot simpler so there is less code maintainance. Most of the work for beammaker is figuring out the formats and making a nice user interface anyways.
import numpy as np
nparticles = 100
mean = [0,0] # mean is zero
cov = C["amplitude"] * C["beta"] * C["muon_mass"] / C["p_z"], 0],<br />[0, C["muon_mass"] * C["amplitude"] / C['beta'] / C['p_z'] # covariance matrix from Penn (thanks Tim!)
x,x_prime = np.random.multivariate_normal(mean,cov,nparticles).T
done. I assume alpha is zero. This uses the super fancy thing I found on wikipedia and had somebody teach me today.
So Tim: do you think this is easier to maintain and/or more clear? If so, feel free to use it. If not, ignore me.
Updated by Rogers, Chris almost 12 years ago
Right - CLHEP has this algorithm and it's what I always used. Hit-miss is bizarre way to do it and probably buggy.
Updated by Carlisle, Timothy almost 12 years ago
- Assignee changed from Rogers, Chris to Carlisle, Timothy
Corrections, after meeting with Chris R:
1) Hit/Miss mc to be removed and replaced with multivariate gaussian.
2) Add basic comments in ConfigDefaults
3) replace pencilbeam datacard with mode, accepting an integer
4) adapt to use Beta_x,y and Emittance_x,y etc as datacards
5) list ConfigDefaults in Birth comments
6) PID as string or integer. E.g. mu- / 13
7) Change variable names according to SpillSchema.py (add Time field)
8) remove Icool format
9) zPosition - fix!!
10) Death returns true
11) comments: matched (multivariate etc), matched how? const Beta, alpha 0 etc
12) Tests: can use xboa to calc. covariance matrix
13) SigmaEnergy, Energy-Time Correlation
Updated by Rogers, Chris almost 12 years ago
How's this going? Would be great to close this issue shortly...
Updated by Carlisle, Timothy almost 12 years ago
Will be done by middle of next week, i.e. before the next meeting.
Updated by Carlisle, Timothy almost 12 years ago
Three running modes are now implemented: matched beta fn, pencil beam & decoupled x,y (can't thing of a decent name!).
remaining issues:
1) test covariance matrices with xboa
2) add basic card tests for mode 2 & 3 (non-zero etc)
3) SigmaEnergy, Energy-Time Correlation
extras:
4) Amplitude momentum correlation
5) dispersion?
6) 6D emittance
Updated by Carlisle, Timothy almost 12 years ago
Ideas:
1) The tests file is getting rather large as more modes are being added, it may make more sense to split it according to each mode i.e. test_MapPyBeamMaker_mode1.py etc - does this sound reasonable? I want to exhaustively test each card for each mode.
2) I've used a dictionary to specify the appropriate cards for each mode:
self.list_of_datacards = {'1' : ['emittance_4D', 'sigma_pz','sigma_time', 'b_z'] ,
'2' : ['sigma_time','sigma_pz'],
'3' :['emittance_x', 'alpha_x', 'beta_x', 'emittance_y',
'alpha_y', 'beta_y', 'sigma_time','sigma_pz'] }
Birth then cross checks this list against the supplied cards. If they match, it is checked (non zero etc). If the user specified it by mistake (e.g emittance_x in mode 1 (matched beta function, use emittance4d card)) then Birth returns false.
I've used a dictionary to try to avoid duplicate code & to streamline it a bit, it probably needs a bit more tidying up though.
3) At the moment all possible cards are initialized and set to zero. My knowledge of python lets me down here...can I simply create them as necessary in Birth? Saves me having superfluous variables.
4) Pylint doesn't like me (I'm heartbroken!), due to lots of branches, lots of variables (addressed in 3) ), and lots of except statements without types (ValueError etc). A sit down with either Chris will prob sort these bits out, down to my python inexperience I imagine.
Updated by Rogers, Chris almost 12 years ago
Are you ready for me to have another look? If so, can you change the assignee to rogers?
Updated by Carlisle, Timothy almost 12 years ago
- Assignee changed from Carlisle, Timothy to Rogers, Chris
Tidying it up atm, found a few schoolboy errors.
Updated by Carlisle, Timothy almost 12 years ago
Done. My connection is pretty slow down here in Devon, but managed to tidy some bits up.
Updated by Rogers, Chris almost 12 years ago
I had a problem checking out your development branch. bzr failed on command like
bzr checkout lp:~t-carlisle1/maus/Devel
with
bzr: ERROR: No such file: ('emrlayer9bar9.dat-20110331150457-pnvhq3en9f9j50xl-5713', 'vassil.verguilov@gmail.com-20110331150634-y5r9larh85u5c2em')
I checked that I could checkout a different (randomly selected) branch. Do you get this error if you attempt to checkout a fresh copy?
Updated by Rogers, Chris almost 12 years ago
- Assignee changed from Rogers, Chris to Carlisle, Timothy
Updated by Rogers, Chris almost 12 years ago
Looks like it - probably some bzr bug, though a bit worrying. I would back up your existing code (just copy to a new directory for safety) then try making a new branch following the instructions here Bzr_usage; then copy any changes across and commit. I wanted to look at your code before making a meeting later in the week...
Updated by Carlisle, Timothy almost 12 years ago
- Assignee changed from Carlisle, Timothy to Rogers, Chris
Ok I've committed it to a new branch: 'test'. Hopefully this one will work...
Updated by Rogers, Chris almost 12 years ago
- Needs pylint running over it
birth()
- Better to do
try:
...
except:
return False
for the whole function? Then you can just return False if you get an exception anywhere.
- assert is redundant. Should get caught in try ... except... statement
- Lots of cut and paste/dead code
key = ["x", "y"]
for i in key:
_code_
key = ["u", "v"]
for i in key:
_code_
Better to do
key = ["u", "v", "x", "y"]
for i in key:
_code_
I would like folks to use the ErrorHandler so we get humane error messages coming through when stuff fails.
try: ... except: return False
for the whole function? Then you can just return False if you get an exception anywhere.
key = ["x", "y"] for i in key: _code_ key = ["u", "v"] for i in key: _code_
Better to do
key = ["u", "v", "x", "y"] for i in key: _code_
I would like folks to use the ErrorHandler so we get humane error messages coming through when stuff fails.
- mode 1, 2, 3 is non-intuitive. I know this is what Mark did, I think it's horrid. Surely better to have mode "pencil", "twiss", "penn"... we aren't FORTRAN coders!
process()
- I like the fact you broke the process function into sub functions. Logic looks correct
tests
- looks okay.
So this is almost ready to go, but it's sort of on the critical path now - so can I be a real pain and do the final polish?
Updated by Rogers, Chris almost 12 years ago
- Target version changed from Future MAUS release to MAUS-v0.1.0
Updated by Rogers, Chris almost 12 years ago
- Assignee changed from Rogers, Chris to Carlisle, Timothy
Okay, I had an iteration - now in my branch. Quite a few changes:
- I moved the stuff that samples particles from parent distribution into a subclass called beam. This means that we can do the next thing:
- I added functionality that allows for multiple beam definitions; this means we can (for example) define separate distributions for electron, pion, muon beams. This wasn't in the original brief but is something we need
- I broke down the structure you had so that the birth() procedure sets up the beam ellipse (covariance matrix). Just seems a more robust way to do it. Then we can just call numpy multivariate gaussian algorithm indepently on some beam matrix
- I added functionality to make a reasonably realistic time distribution, so sawtooth or uniform. This wasn't in the original brief but is something we need
- I extended the test coverage, documentation.
So it's in
https://code.launchpad.net/~chris-rogers/maus/devel
available by doing
bzr branch lp:~chris-rogers/maus/devel
The following source files
src/map/MapPyBeamMaker/*.py
are used
Tim, can you look at the code/do a review? I might also ask Matt Littlefield to have a look. Would like to push it by end of next week.
Cheers!
Updated by Rogers, Chris almost 12 years ago
Nb: I notice there is a test fail in my branch for some stuff I'm working on for the ErrorHandler. Sorry about that, it doesn't affect any of the code here.
Updated by Littlefield, Matthew almost 12 years ago
I've looked over the code and although I don't really understand the nitty gritty of the physics involved I get the jist and it looks ok. The code itself is good maybe a little untidy in places where lines need to be split so you don't have long lines of code but it can still be followed. :-)
Updated by Tunnell, Christopher almost 12 years ago
Where is it untidy? Providing line numbers means that we can improve constructively.
For the physics that you don't understand (which I imagine is a particle/accelerator divide) can Chris provide in the comments a book recommendation? Maybe a section in Ted Wilson's book or something?
Updated by Rogers, Chris almost 12 years ago
I realise I forgot to fill out some of the equations at the top of beam.py
Some references is a good idea
Updated by Rogers, Chris over 11 years ago
- Status changed from Open to Closed
- % Done changed from 0 to 100
Committed in r634
Updated by Rogers, Chris over 11 years ago
- Target version changed from MAUS-v0.1.0 to MAUS-v0.0.7