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 ...";
music_message.flush("warning");
}
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;
music_message.flush("warning");
}
// 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.";
music_message.flush("warning");
} 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...";
music_message.flush("info");
stringstream tmpStr;
tmpStr << "./outputs/vnchdeta_pT_" << pT_min << "_" << pT_max << ".dat";
ofstream outfile;
outfile.open(tmpStr.str().c_str(),ios::trunc);
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;
}
outfile.close();
}
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
Warning
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.