Project

General

Profile

Feature #454

Review beam maker code

Added by Rogers, Chris over 12 years ago. Updated about 12 years ago.

Status:
Closed
Priority:
Normal
Category:
Simulation
Target version:
Start date:
17 May 2011
Due date:
30 May 2011
% Done:

100%

Estimated time:
Workflow:

Description

Just want to remind myself to have a look


Related issues

Related to MAUS - Feature #357: Beam makerClosedCarlisle, Timothy01 March 2011

Actions
#1

Updated by Rogers, Chris over 12 years ago

Relevant files:

MapPyBeamMaker
test_MapPyBeamMaker

MapPyCreateSpill
test_MapPyCreateSpill

#2

Updated by Rogers, Chris over 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

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
#3

Updated by Rogers, Chris over 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.

#4

Updated by Carlisle, Timothy over 12 years ago

Sure thing.

I rang earlier, but will try again after 2.30pm.

#5

Updated by Rogers, Chris over 12 years ago

Sorry was in London today...

#6

Updated by Carlisle, Timothy over 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

#7

Updated by Tunnell, Christopher over 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.

#8

Updated by Rogers, Chris over 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

#9

Updated by Rogers, Chris over 12 years ago

What is the name of your branch? So I can check it out

#10

Updated by Tunnell, Christopher over 12 years ago

bzr info

#11

Updated by Rogers, Chris over 12 years ago

(Tim needs to do it though)

#12

Updated by Tunnell, Christopher over 12 years ago

said it for him so he knew how to complete the request ;)

#13

Updated by Tunnell, Christopher over 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.

#14

Updated by Rogers, Chris over 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.

#15

Updated by Carlisle, Timothy over 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

#16

Updated by Rogers, Chris over 12 years ago

How's this going? Would be great to close this issue shortly...

#17

Updated by Carlisle, Timothy over 12 years ago

Will be done by middle of next week, i.e. before the next meeting.

#18

Updated by Carlisle, Timothy over 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

#19

Updated by Carlisle, Timothy over 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.

#20

Updated by Rogers, Chris over 12 years ago

Are you ready for me to have another look? If so, can you change the assignee to rogers?

#21

Updated by Carlisle, Timothy over 12 years ago

  • Assignee changed from Carlisle, Timothy to Rogers, Chris

Tidying it up atm, found a few schoolboy errors.

#22

Updated by Carlisle, Timothy over 12 years ago

Done. My connection is pretty slow down here in Devon, but managed to tidy some bits up.

#23

Updated by Rogers, Chris over 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?

#24

Updated by Rogers, Chris over 12 years ago

  • Assignee changed from Rogers, Chris to Carlisle, Timothy
#25

Updated by Carlisle, Timothy over 12 years ago

same error...is my branch at fault?

#26

Updated by Rogers, Chris over 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...

#27

Updated by Carlisle, Timothy over 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...

#28

Updated by Rogers, Chris over 12 years ago

Few comments:
  • 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.

  • 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?

#29

Updated by Carlisle, Timothy over 12 years ago

yep all yours!

#30

Updated by Rogers, Chris over 12 years ago

  • Target version changed from Future MAUS release to MAUS-v0.1.0
#31

Updated by Rogers, Chris about 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!

#32

Updated by Rogers, Chris about 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.

#33

Updated by Littlefield, Matthew about 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. :-)

#34

Updated by Tunnell, Christopher about 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?

#35

Updated by Rogers, Chris about 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

#36

Updated by Rogers, Chris about 12 years ago

  • Status changed from Open to Closed
  • % Done changed from 0 to 100

Committed in r634

#37

Updated by Rogers, Chris about 12 years ago

  • Target version changed from MAUS-v0.1.0 to MAUS-v0.0.7
#38

Updated by Rogers, Chris about 12 years ago

  • Tracker changed from Support to Feature

Also available in: Atom PDF