Contents

Step 2: Generate data samples

Warning

The pycompwa is no longer maintained. Use the ComPWA packages QRules, AmpForm, and TensorWaves instead!

In this section, we will use the amplitude model that we created with the expert system to generate a data sample via hit & miss Monte Carlo. This requires us to work with the User Interface to ComPWA, the C++ back-end of pycompwa.

2.1 Generate

In the previous step, we created an amplitude module and exported it as an XML file. We now need to go through a few steps to generate a phase space sample and a data sample for the decay \(J/\psi \rightarrow \gamma\pi^0\pi^0\):

  1. Import the section of particle definitions that is embedded in the XML file.

  2. Build a Kinematics object following the specifications of the dynamics in the XML model file (in this case, the helicity formalism). The second argument of this function is a particle list.

  3. The information in the Kinematics instance is sufficient for generating a phase space sample. This sample defines the space on which to evaluate the intensities of the amplitude model.

  4. For generating a data sample, you require an Intensity profile for this specific decay. You can construct such an object from the model XML file using the IntensityBuilderXML. Because we want the intensities to be normalized, you have to pass the generated phase space sample as an argument to the intensity builder. Note that during the creation of the Intensity object, the Kinematics instance is updated with the subsystems required for the decay topology.

Now all building blocks for generating our data sample are at hand and you can define a data sample of arbitrary size!

import pycompwa.ui as pwa  # interface to ComPWA

2021-10-27 10:27:26,316 INFO [default] Logging to file disabled!
2021-10-27 10:27:26,316 [INFO] Log level: INFO
2021-10-27 10:27:26,316 [INFO] Current date and time: Wed Oct 27 10:27:26 2021

# specify kinematics
particle_list = pwa.read_particles("model.xml")
kinematics = pwa.create_helicity_kinematics("model.xml", particle_list)

# specify generators
generator = pwa.EvtGenGenerator(
    kinematics.get_particle_state_transition_kinematics_info()
)
random_generator = pwa.StdUniformRealGenerator(12345)

# generate phase space sample
phsp_sample = pwa.generate_phsp(int(1e5), generator, random_generator)

# create intensity profile
intensity_builder = pwa.IntensityBuilderXML(
    "model.xml", particle_list, kinematics, phsp_sample
)
intensity = intensity_builder.create_intensity()

# generate data
data_sample = pwa.generate(
    int(1e4), kinematics, generator, intensity, random_generator
)

2021-10-27 10:27:26,709 [INFO] HelicityKinematics::HelicityKinematics() | Initialized kinematics for reaction ( J/psi )->( gamma[ID=2] pi0[ID=3] pi0[ID=4] )
Event position to final state ID mapping:
0: 2
1: 3
2: 4

2021-10-27 10:27:26,711 [INFO] Generating phase-space MC: [100000 events] 

?25l 
 0.00min [ ---------- 0.00% ] 0.00min , end time: 10:27AM               ?25h

?25l 
 0.00min [ ********** 100.00% ] 0.00min , end time: 10:27AM               ?25h

2021-10-27 10:27:26,794 [INFO] Setting phase space sample weights...

2021-10-27 10:27:26,802 [INFO] Updating data container content...

2021-10-27 10:27:28,496 [INFO] Generating hit-and-miss sample: [10000 events] 

?25l 
 0.00min [ ---------- 0.00% ] 0.00min , end time: 10:27AM               ?25h

2021-10-27 10:27:29,082 [INFO] Tools::generate() | Error in HitMiss procedure: Maximum value of random number generation smaller then amplitude maximum! We raise the maximum to 9.92631 value and restart generation!

?25l 
 0.00min [ ---------- 0.00% ] 0.00min , end time: 10:27AM               ?25h

2021-10-27 10:27:31,090 [INFO] Successfully generated 10000 with an efficiency of 0.0740741
?25l 
 0.03min [ ********** 100.00% ] 0.00min , end time: 10:27AM               ?25h?25l 
 0.03min [ ********** 100.00% ] 0.00min , end time: 10:27AM               ?25h

Note

pycompwa.ui is the python interface to ComPWA’s C++ modules. It has four major important components for evaluating intensities:

  • The EventCollection class is simply an event-based list of four-momenta with a list of Particle IDs.

  • The Kinematics class allows one to compute kinematic variables, like \((s,\theta,\phi)\) in the helicity formalism, from the momentum tuples in an EventCollection.

  • The DataSet class is the result when you convert() an EventCollection. It is essentially a table of kinematic variables for each event in that original collection of four-momenta.

  • The Intensity class allows one to compute an intensity (an absolute value) for each data point in the DataSet. You do this computation with the evaluate() method (see Step 4: Analyze results).

Hint

For more complicated data sample structures, see Importance sampling.

2.2 Exporting and importing

The pycompwa.ui module allows you to export the data samples to two different file formats: ROOT files and ASCII files.

# to ROOT
pwa.write_root_data(data_sample, output_file="generated_data.root")
pwa.write_root_data(phsp_sample, output_file="generated_phsp.root")
# to ASCII
pwa.write_ascii_data(data_sample, output_file="generated_data.dat")
pwa.write_ascii_data(phsp_sample, output_file="generated_phsp.dat")

2021-10-27 10:27:31,265 [INFO] Root::writeData(): writing vector of 10000 events to file generated_data.root

2021-10-27 10:27:31,405 [INFO] Root::writeData(): writing vector of 100000 events to file generated_phsp.root

2021-10-27 10:27:31,999 [INFO] Logging to file disabled!

The syntax for importing parallels that:

imported_data = pwa.read_root_data(input_file="generated_data.root")
imported_phsp = pwa.read_ascii_data(input_file="generated_phsp.dat")

print("Imported data events:", len(imported_data.events))
print("Imported phsp events:", len(imported_phsp.events))

first_event = imported_data.events[0]
print("First data event:")
print(first_event.four_momenta)
print(first_event.weight)

Imported data events: 10000
Imported phsp events: 100000
First data event:
FourMomentumList[(-0.520903,0.885259,1.21872), (0.653672,-0.813022,1.46359), (-0.132769,-0.0722372,0.414596)]
1.0

Note that the ASCII (.dat) files contain a header section that informs you about the final state of the file. You will have to prepend this section to the file yourself if you want to use an ASCII file that was exported from another framework.

with open("generated_data.dat") as f:
    for _ in range(5):
        print(f.readline()[:-1])

<header>
	Pids: [22, 111, 111]
	Order: px py pz E
	Unit: GeV
</header>