5. Final hadron dynamics: decays & afterburners

Particlization describes the conversion of the fluid to hadronic degrees of freedom. In practice, in MUSIC, this means that the momentum distribution for each species of hadrons is tabulated from the particlization hypersurface.

Physically particlization is the transition from the strongly-coupled quark-gluon plasma reconfined to a gas of hadrons. These hadrons continue to interact until it is too dilute. Moreover most of the hadrons in this gas are unstable and decay into stable and unstable hadrons. These decays continue until only hadrons that are stable [*] remain.

[*]In the present case, stability means that their lifetime is sufficient to reach the detectors.

MUSIC does not include a simulation of hadronic interaction. Including such interaction requires the addition of hadronic transport models such as UrQMD and SMASH.

MUSIC does include the decay of unstable hadrons however. In general these decays have a much larger effect of the final momentum distribution of hadrons that hadronic interactions.

To compute hadronic decays, MUSIC reads the output of mode 3, the “yptphiSpectra.dat”, and read the branching ratio and decay channels from the “pdg05.dat” or corresponding file discussed in the previous section (see box ” Species in Cooper-Frye”). The code performing the decay can be found in source file “reso_decay.cpp”

MUSIC produces a new file, “FyptphiSpectra.dat”, which has exactly the same format as the corresponding file in Cooper-Frye: a tabulation of the momentum distribution of each species of hadrons considered.

5.1. Comparing with data: calculating observables

After hadronic decays are calculated, the tabulated momentum spectra in “FyptphiSpectra.dat” can be compared with experimental data. The last step is to compute observables.

A selection of observables is already calculated in MUSIC within function “void Freeze::compute_final_particle_spectra_and_vn([…])” in source file “freeze_pseudo.cpp”.

void Freeze::compute_final_particle_spectra_and_vn(InitData* DATA) {
    // take tabulated post-decay spectra and
    // compute various observables and integrated quantities
    int status = mkdir("./outputs", 0755);   // output results folder
    if (status != 0) {
        music_message << "Can not create folder ./outputs ...";
    ReadSpectra_pseudo(DATA, 1, 1);  // read in particle spectra information

    if (DATA->NumberOfParticlesToInclude > 200) {
        if ((DATA-> whichEOS != 7 && particleMax < 320)
                || ((DATA->whichEOS == 7) && particleMax < 327)) {
            music_message << "Not all hadronic species are computed "
                          << "particleMax = " << particleMax;
        // calculate spectra and vn for charged hadrons
        Output_charged_hadrons_eta_differential_spectra(DATA, 1, 0.01, 3.0);
        Output_charged_hadrons_pT_differential_spectra(DATA, 1, -0.5, 0.5);
        Output_charged_IntegratedFlow(DATA, 0.01, 3.0, -0.5, 0.5);
        // for PHENIX pT cut
        Output_charged_IntegratedFlow(DATA, 0.15, 3.0, -0.5, 0.5);
        // for ALICE or STAR pT cut
        Output_charged_IntegratedFlow(DATA, 0.2, 3.0, -0.5, 0.5);
        // for CMS pT cut
        Output_charged_IntegratedFlow(DATA, 0.3, 3.0, -0.5, 0.5);

        Output_charged_IntegratedFlow(DATA, 0.4, 3.0, -0.5, 0.5);
        // for ATLAS pT cut
        Output_charged_IntegratedFlow(DATA, 0.5, 3.0, -0.5, 0.5);

    if (particleMax < 21) {
        music_message << __func__
                      << "cannot compute pion, kaon and proton observables, "
                      << "some of the spectra are not available.";
    } else {
        // calculate spectra and vn for identified particles
        int pid_list [] = {211, -211, 321, -321, 2212, -2212};
        int pid_list_length = sizeof(pid_list)/sizeof(int);
        for (int k = 0; k < pid_list_length; k++) {
            int number = pid_list[k];
            int id = partid[MHALF+number];
            if (id < particleMax) {
                OutputDifferentialFlowAtMidrapidity(DATA, number, 1);
                OutputDifferentialFlowNearMidrapidity(DATA, number, 1);
                OutputIntegratedFlow_vs_y(DATA, number, 1, 0.01, 3.0);

As you can see above, function “compute_final_particle_spectra_and_vn()” actually called other functions that compute the observables. An example of function where observables are calculated is:

// Output pT and phi integrated charged hadron spectra as a function of eta
void Freeze::Output_charged_hadrons_eta_differential_spectra(
        InitData *DATA, int full, double pT_min, double pT_max) {
    music_message << "Collecting charged hadrons dN/deta and vn vs eta...";
    stringstream tmpStr;
    tmpStr << "./outputs/vnchdeta_pT_" << pT_min << "_" << pT_max << ".dat";

    ofstream outfile;

    outfile << "#eta  dNch/deta  v1cos  v1sin  v2cos  v2sin  v3cos  v3sin  "
            << "v4cos  v4sin  v5cos  v5sin  v6cos  v6sin  v7cos  v7sin\n";

    double eta_min = DATA->dNdy_eta_min;
    double eta_max = DATA->dNdy_eta_max;
    int nrap = DATA->dNdy_nrap;
    double deta = (eta_max - eta_min)/(nrap - 1);
    int max_flow_order = 7;
    for (int ieta = 0; ieta < nrap; ieta++) {
        double eta_local = eta_min + ieta*deta;
        double dNch = get_Nch(DATA, pT_min, pT_max, 0, eta_local, eta_local);
        outfile << eta_local << "  "  << dNch << "  ";
        for (int iorder = 1; iorder < max_flow_order; iorder++) {
            double *vn_temp = new double [2];
            get_vn_ch(DATA, pT_min, pT_max, 0, eta_local, eta_local,
                      iorder, vn_temp);
            outfile << vn_temp[0] << "  " << vn_temp[1] << "  ";
            delete [] vn_temp;
        outfile << endl;

Function “Output_charged_hadrons_eta_differential_spectra()” itself calls helper functions like “get_Nch()” and “get_vn_ch()” that return observables like the multiplicity and \(v_n\) of charged hadrons. The results is outputted in a text file saved in the directory “./outputs/”.

The input of the helper functions like “get_Nch()” is the following:

// Return charged hadron yield in specified range of phase space
// You'll almost certainly want to set yflag=0 to calculate
// over a fixed range in eta
// if minrap==maxrap, returns dN/eta.  If minpt==maxpt, dN/dpt.
// If both, dN/dpt/deta.  Otherwise, total yield N
double Freeze::get_Nch(InitData *DATA, double minpt, double maxpt, int yflag,
                       double minrap, double maxrap) { [...]  }

An observable like the number of charged hadrons is integrated over \(\phi\). This leaves two momentum degrees of freedom, \(p_T\) and either rapidity \(y\) or pseudorapidity \(\eta\). The number of charged hadrons can be integrated over a range of \(p_T\) and \(y/\eta\), or be given at a specific value of \(p_T\) and \(y/\eta\). As explained in the function comments, “get_Nch()” can do both, depending on its parameters. The parameters of “get_Nch()” are:

  • “minpt” & “maxpt”: either the range of \(p_T\) to integrate over, or the value of \(p_T\) are which \(dN/dp_T\) is requested
  • “yflag”: whether rapidity \(y\) or pseudorapidity \(\eta\) should be used
  • “minrap” & “maxrap”: either the range of :\(y/\eta\) to integrate over, or the value of \(y/\eta\) are which \(dN/dy\) or \(dN/d\eta\) is requested


Since these hadronic observables are computed from a pre-tabulated spectra, “minpt/maxpt/minrap/maxrap” should all be within the range of the tabulated momentum distribution (i.e. content of “FyptphiSpectra.dat”, which uses the same grid as Cooper-Frye, “yptphiSpectra.dat”)

Similarly, since the momentum is discretized with pseudo-rapidity, any value of rapidity requested should be calculable from the available discretized table.

Computing most observables by modifying the existing code and by reference to the above examples should be straightforward.