info about whizard commands

Asked by Andrea Paolo Puppin

Hi ,

sorry to bother you.

I’m a student of physic at the university of study in Milan.

i wanna simulate two beams of electrons and photons to generate two electron and one positron( triplet pair production).

I started to use whizard but i have some problems to understand some of the output.

I wanna kindly ask you( i wanna simulate two colliding beams ) :

1)is it possible to set the momentum spread of my two beams?

2)is it possible to set some parameter for considering the physical size of the two beams ?

3) can i ( in only one simulation) consider a lot of angles between the incident particles of the two beams or i have to do a single simulation for every angle(and then eventually merge all the output for every single simulation)

4)the observable which whizard produce are in the centre of the laboratory or in the centre of mass?

5)what’s the command to produce an output file(it should be the analysis.dat file) where i put in every colomn all the observable for every single particles?

Thanks you in advance and sorry for all these questions,

Kind regards

Question information

Language:
English Edit question
Status:
Solved
For:
WHIZARD Edit question
Assignee:
Juergen Reuter Edit question
Solved by:
Juergen Reuter
Solved:
Last query:
Last reply:
Revision history for this message
Juergen Reuter (j.r.reuter) said :
#1

Hi Andrea,
thanks for reaching out to us. Launchpad is always better than via mailing list, and you found it yourself before we could answer.
1)
It is possible to use a Gaussian beam spread on the energies of the two beams, using
beams = e1, E1 => gaussian
gaussian_spread1 = 0.01%
gaussian_spread2 = 0.01%
with corresponding spreads of the two beams, respectively. At the moment, a spread of the momenta or the angles themselves is not foreseen. This was already asked for in the past once, but this has not been turned into an official feature request. We will discuss it at our next collaboration meeting.
2)
No, physical sizes of beams or impact parameters are not foreseen. This has to be done via a dedicated beam simulation, e.g. via Guinea Pig.
3)
Whizard does a dedicated phase space generation for every kinematic configuration, and a different beam angle constitutes a different kinematic configuration. You can, however, semi-automatize this by doing a scan, along the lines of something like
beams = e1, E1
beams_momentum = 50 GeV, 50 GeV
n_events = 1
sample_format = lhef
scan real x (2 degree => 11 degree /+ 1.5 degree) {
  beams_theta = 0, x
  integrate (my_process)
 $sample = sprintf "eeff_%1.3f" (real (x))
  simulate (my_process)
}
This puts everything into different event files for each crossing angle. Then you can use tools to merge event samples (e.g. Yoda or Root). There is, in principle, the possibility to use the alt_setup feature of Whizard to reweight results of one integration to another one with different parameters. As the phase space integration is left unchanged, I would not use it here because it might lead to systematic errors, as the phase space and the Lorentz system changes with crossing angle. (A dedicated beam simulation tool like Guinea Pig has only a few standard candle physics processes with closed formulae for there cross sections, so it might be able to easier scan over beam-specific parameters changing the Lorentz frame.)
4)
Whizard produces results in the laboratory frame, which is assumed to be the one of the collider.
5)
It is not clear to me exactly what you want to do. Do you want to generate several histograms at the same time, or calculate observable over an event sample?
Please give us a bit more information, especially on the last point.
Cheers,
    JRR (Juergen)

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#2

hi,
first of all thank you for your answers.
For the question number 5 i mean:
5)
how can i generate an output file where for every row i have a columns with all observable of the particle of my interest(like angle, energy,momentum... ) because i wanna see for every positon all his observables.
what is also not clear to me is this:
6)i'm trying to study triplet pair production(so my outputs particles are e1,E1,E1), how can i do what i asked in question 5 for the two electrons, how can whizard recognize one electron from the other?
7)what exactly mean if i set a gaussian beam? i mean, whizard does the integration considering a single process at a time(in my problem it consider a colliding photon and electron) and it calculates the cross section for n events.
the problem is that i wanna study triplet pair production for two colliding beams but if i understood well whizard can't exactly simulate this.
i must do multiple whizard run and then eventually merge in some way, correct?
Thanks again for you time,
Kind regards

Revision history for this message
Juergen Reuter (j.r.reuter) said :
#3

Your question 5) (and 6)) sounds like you are looking for a specific event format. E.g. LHE or HepMC are ASCII files that list momenta, energies, masses and helicities of particles in the event output. If you have different particle-wise observables this sounds again like something that best a reconstruction or analysis software tool delivers.
In a Whizard analysis you can get a specific particle observable by defining e.g.
analysis = record pt_e1 (eval Theta [extract index 1 [e1]]]);
                    record pt_e2 (eval Theta [extract index 2 [e1]]])
This would give you really the angle of the first and second electron in the process as ordered by the MC (which is of course completely unphysical). Usually this is used on an ordered list to given the hardest, 2nd-hardest, etc. lepton or jet.
You can see an example in share/examples/LHC_VBS_likesign.sin.
7) A Gaussian beam is one where the energies of the incoming beam particles are not monochromatic, but smeared by a Gaussian with a standard deviation given usually as percentage of the mean of the beam energy.

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#4
Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#5

Thanks again, but unfortunately i have other questions.
1) i'm trying to run this code but in the analysis.dat file is empty.
process ee = e1, A => e1,e1, E1

beams = e1, A => gaussian

gaussian_spread1= 0.1%

gaussian_spread2= 0.2%

sqrts= 10 MeV

write_analysis
integrate (ee)

analysis= record pt_e1(eval Theta[extract index 1 [e1]]) ;

          record pt_e1(eval Theta[extract index 2 [e2]])

simulate (ee) {

n_events=1000

 }

 compile_analysis

also, i havent' well understood:

2) whizard create by default a file where i can find all observables?
3) i have tried to create a file analysis.dat with the code written below.
this create an analysis.dat file with the first 10000 raws with all Px,Py of the photons( using plots which is probably wrong).
The following 10000 raws are with All E,Theta always of the photons.
what i wanted to ask in previous question is if i can have only 10000 raws with all the quantities i wanna study,
process eeprovaluca = A, e1 => e1,A

beams = A, e1=> gaussian

gaussian_spread1= 25.0%

gaussian_spread2= 0%

beams_momentum = 10 MeV , 10 MeV

integrate (eeprovaluca)

plot Px_Py

plot pippo

analysis= record Px_Py (eval Px["A"], eval Py["A"]) and

   record pippo (eval E["A"], eval Theta["A"])

!record cose1 (eval cos (Theta) ["e1"])

!analysis= record prova (eval Px["E1"], eval Px["e1"])

write_analysis

simulate (eeprovaluca) {

n_events=10000

 }

 compile_analysis

Revision history for this message
Juergen Reuter (j.r.reuter) said :
#6

1)
Sorry, in my example I forgot to define the observables, you have to do something like
histogram pt_e1 (0,4,1)
histogram pt_e2 (0,4,1)
analysis= record pt_e1(eval Theta[extract index 1 [e1]]) ;
   record pt_e2(eval Theta[extract index 2 [e1]])

It might also be useful for you to study the manual, especially Chap. 4 and 5.
2)
Whizard by default writes out an event files with all the particle objects and their kinematic information. This is, however, a machine-specific file, with ending .evx. You need to ask for an event file with format LHE (ending .lhe) e.g., by setting sample_format = lhef (or, alternatively, if linked HepMC, ending .hepmc. Those contain the particles in ASCII format.
3)
You can actually write out up to 4 observables into 4 different columns, using something like

plot Px_Py
analysis= record Px_Py (eval Px["A"], eval Py["A"], eval E["A"], eval Theta["A"])

The corresponding PDF would not be correct, but you can just use
write_analysis { out_file = "foo.dat" }
to only create the .dat file.

Note that the internal Whizard analysis is not very sophisticated, and it might be beneficial for you to use something like Rivet or Root for an analysis.
You can also several .dat files and merge their data with corresponding tools.

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#7

okay thank you very much !!!
i'll probably write again when i'll try to do the simulation for a lot of angles(i'll try the code you wrote )
i just wanna ask you one more thing
1) i tried to read the manual but what i didn't understood is how exactly whizard works, i mean how (not all the theory but just the idea) it calculates the cross section of the progress.
when i set nevents=1000 (for like e1 + A = e1 + e1 + positron, so triplet pair production) the program simulate more than 1000 events(considering all the others process during interation ) and not all these events produces the process i wanna study.
so the program generates for exemple more than 1000 events (like 3000 but i think it depence on the cross section of the event i wanna study) until it arrives to generate 1000 triple pair production(in my case) and then it calculates the cross section dividing the total events/nevents?
sorry because it's probably a stupid question,
Kind regards

Revision history for this message
Juergen Reuter (j.r.reuter) said :
#8

Dear Andrea,
there are several aspects to answer your question: first of all, the numbers that you first see appearing on screen is the adaptation of phase space, where Whizard's internal multi-channel Monte Carlo integrator VAMP tries to map out the singularities of phase space as close as possible to a constant function, such that a typical (hit-and-miss-like) MC veto procedure is as efficient as possible. This is called importance sampling. This mapping/adaptation of phase space works the better the more phase space points you try out. These are the numbers that you see, which are phase space calls, but not yet events. (This is like in deep learning where you would take much larger training samples than your expected event numbers). The adaptation of phase space is then encoded in grids which are read in for the event generation. Event generation can be weighted or unweighted, where weighted events are sampled from a uniform random number distribution and give back phase space points with a weight given by the product of the phase space flux factor and the squared matrix element. These are then unweighted, again by a Monte Carlo veto procedure, such that the events finally generated have unit weight but are generated from a random number distribution that are concentrated about areas in phase space with (originally) highest weights. The better the phase space adaptation has worked, the higher is the unweighting efficiency, which is by definition the number of weighted events per single unweighted event to be generated. For a 2->2 process this is typically 10-30%, meaning you need 3-10 weighted events to generate one unweighted event.
Finally, if you specify a process with exclusive final state, there are also irreducible background processes.
Cheers,
   JRR

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#9

process eeprovamultifasci = A, e1 => e1,e1,E1
beams = A, e1 => gaussian
beams_momentum = 50 keV, 10 GeV
n_events = 100
sample_format = lhef
scan real x (2 degree => 11 degree /+ 1.5 degree) {
  beams_theta = 0, x
  integrate (eeprovamultifasci)
 $sample = sprintf "eeff_%1.3f" (real (x))
  simulate (eeprovamultifasci)
}

okay , i hope this one of my last question(sorry)
I tried this code you sent me in your first answers( for triplet pair production study process) but it doesn't work.
what i wanna obtain is a programm which simulate for more angles and it produces more analysis.dat file (where i wanna also put my observable with plot) everyone for a specific angle and then i'll merge them with other tools.
is it possible?
thanks again,
Kind regards

Revision history for this message
Juergen Reuter (j.r.reuter) said :
#10

Sorry, there was a typo in my example, a missing equal sign, it must read:

scan real x = (2 degree => 11 degree /+ 1.5 degree) {

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#11

okay the code works!!
but is it possible to write a code to create a file analysis.dat for every run (maybe colling analysistheta1.dat, analysistheta2.dat....).
or a code to create a single analysi.dat file where (for every run/angle i can plot all observable i wanna study)?
i was trying this :
process eeprovamultifasci = A, e1 => e1,e1,E1
beams = A, e1 => gaussian
beams_momentum = 50 keV, 10 GeV

sample_format = lhef
scan real x= (2 degree => 11 degree /+ 1.5 degree) {
plot angle
plot Px_Py
n_events = 10
  beams_theta = 0, x
  integrate (eeprovamultifasci)
 $sample = sprintf "eeff_%1.3f" (real (x))
 analysis= record Px_Py (eval Px["A"], eval Py["A"]) and
      record angle (eval E["A"], eval Theta["A"])
      write_analysis
  simulate (eeprovamultifasci)
}
compile_analysis
with this code i was trying to create a file analysis.dat where for every run(so for different angles) i have 4 colomuns with px,py,e,theta , i know it becomes a long file but i don't know if there's a sneaky way to do it.
kind regards

Revision history for this message
Best Juergen Reuter (j.r.reuter) said :
#12

I think you need to histogram (then the entries in the data file are filled), or plot single variables as is shown in the example files. They are in the folder share/examples and also documented in the manual. Besides this, I told you that Whizard's internal analysis is not very sophisticated. All the things you mention you could easily encode in an analysis tool like Rivet. I would suggest to close this ticket as 'solved', and in case you have a more specific new question, open a new one. Thank you.

Revision history for this message
Andrea Paolo Puppin (2andre3) said :
#13

Thanks Juergen Reuter, that solved my question.