Relevant files:

MapPyBeamMaker
test_MapPyBeamMaker

MapPyCreateSpill
test_MapPyCreateSpill

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.

Sure thing.

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

Sorry was in London today...

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

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.

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

bzr info

(Tim needs to do it though)

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

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.

Right - CLHEP has this algorithm and it's what I always used. Hit-miss is bizarre way to do it and probably buggy.

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

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

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

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.

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

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

I had a problem checking out your development branch. bzr failed on command like

bzr checkout lp:~t-carlisle1/maus/Devel

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?

same error...is my branch at fault?

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...

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

yep all yours!

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.

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. :-)

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?

I realise I forgot to fill out some of the equations at the top of beam.py

Some references is a good idea