Foreword

The two more common methods for asymmetries extraction are the so called "1D" (1-Dimension) and "UL" (Unbinned Likelihood). Nowadays the usual strategy is the following: the UL is used for the final results, and 1D for the systematics or to do some "quick" tests, due to its simplicity. From the technical point of view, the 1D method is easy to implement, while UL requires more time and a more complex code. This is also translating in different time of execution to have an idea of the order of magnitude, while 1D takes some minutes, the UL takes several hours or days ( depending on which minimizer is used).

Said that, an introduction to the methods can be found in some doctorate theses. E.g. at pags 64-74 of Phd G.Pesaro there is a description of the two methods ( and also of a third one "2D", that can be considered a method in between them). At pag. 59, par. 3.5, there is general description of the principle of the measurement, that is how data are combined.

Unbinned Likelihood

In order to deepen the knowledge about the UL, there is a detailed COMPASS note: 2009-13.

Guidelines for the implementation of the likelihood using MINUIT in ROOT

In the following, some guidelines to implement the UL using MINUIT in the ROOT framework are given. This is one of the recipe (F.Sozzi, Trieste) used for the releases. Other recipes with slighlty different functions (e.g. acceptance terms in) or different minimizers have also been used by other people along the years.

Definition of the likelihood

• Definition of PDF functions for events coming from "spin up" and "spin down" :
  f_pos = 1.+U1 cos(φ_h)+U2 cos(2. φ_h)+A1 sin(φ_h+φ_S-π)+A2 sin(3 φ_h-φ_S)+A3 sin(φ_h-φ_S)+A4 cos(φ_h-φ_S)+A5 sin(φ_S)+A6 sin(2. φ_h-φ_S)+A7 cos(φ_S)+A8 cos(2 φ_h-φ_S);
  f_neg = 1.+U1 cos(φ_h)+U2 cos(2. φ_h) - A1 sin(φ_h+φ_S-π) - A2 sin(3 φ_h-φ_S) - A3 sin(φ_h-φ_S) - A4 cos(φ_h-φ_S) - A5 sin(φ_S) - A6 sin(2. φ_h-φ_S) - A7 cos(φ_S) - A8 cos(2 φ_h-φ_S));
  Note that: A1 = A_Col; A3 = A_Siv. The π term in the Collins angle is due to the COMPASS sign convention, that is opposite to the Trento convention.
  The free parameters are 10: [U1-U2] are in the unpolarized part, the eight [A1-A8] are the transverse spin asymmetries, that change sign with the spin orientation.

• Definition of a PDF function for each cell.
  In this example, we are using 4 cells like in the 1D double ratio case. Note that you can define 8 cells as well, like in "quadrupole ratio" case.
  So, naming the two periods with opposite polarisations as p1 and p2, the PDFs are:
  f_pos_p1 = c0 * f_pos
  f_pos_p2 = c1 * f_pos
  f_neg_p1 = c2 * f_neg
  f_neg_p2 = c3 * f_neg
  where [c0,c3] are 4 additional free parameters.
  In total the free parameters are 10 + 4 = 14 (or 18, in case of "QR" implementation)

• Creation of the function to minimize in MINUIT (FCN).
  Instead of maximising the likelihood, one minimizes its logarithm, multiplying it by -1.
  Within a loop on the hadrons, the 4 logarithms corresponding to the 4 cells are built:
  lnLpos_p1 += log(f_pos_p1);
  lnLpos_p2 += log(f_pos_p2);
  lnLneg_p1 += log(f_neg_p1);
  lnLneg_p2 += log(f_neg_p2);
  
  After the loop on the hadrons, the FCN function can be built:
  
  FCN = - lnLpos_p1/n1_pos - lnLneg_p1/n1_neg - (other 2 pieces for 2nd period) +
  + Int_f_pos_p1/n1_pos + Int_f_neg_p1/n1_neg + (other 2 pieces for 2nd period)
  
  where n1_pos and n1_neg are the number of hadrons with up and down polarisation in period p1.
  The parts Int_f_pos_p1, Int_f_neg_p1 are the integrals over [-π , π] in φ_h and φ_S of the functions f_pos_p1 and f_neg_p1 respectively
  This is the part of "normalization" of the likelihood, indicated in the standard UL formula as the "exp (-N)" term that multiplies each of the 4 terms in the likelihood (see formula 3.18 in Giulia Pesaro's thesis).
  In case of perfect acceptance without "holes" in phih and phis, the integrals of the above functions are then simply:
  Int_f_pos_p1 = c0 * 4*π²
  Int_f_pos_p2 = c1 * 4*π²
  Int_f_neg_p1 = c2 * 4*π²
  Int_f_neg_p2 = c3 * 4*π²
  For a first implementation of the code, one can use this formula and check if it is working.
  For the final implementation of the code, one should consider the actual acceptance integrating only over the range in phih and phis in which there is statistics. In any case note that the effect in "standard" bins and configuration is small, since usually there are no holes in our acceptance.
  Here is an example of how the integral can be evaluated numerically. Do a loop on phih and phis values, dividing the range [-π , π] in pieces (e.g. 100 times 100 "bins". Note that this division is arbitrary) .
  Evaluate the PDFs using the mean value of phih and phis in that bin and multiply for the size of the bin (2pi/100)*(2pi/100).
  This is the contribution of that bin to the integral. To have the total integral over [-π , π] , all contributions should be summed up (one should find again the cx * 4*π² values). To have the integral in a sub-range of the phih and phis space, only the corresponding bins should be added.
  In particular, in this case we are interested to reject the empty bins, in which there is no statistics. In order to know which bins are empty, one can fill at the beginning at the program (before the minimization) a 2D plot in phih and phis with the same binning.

Details about the Minimization procedure.

The MINUIT package allows to personalise the minimization procedure and the algorithms to use for this purpose and for the error evaluation. Here is one recipe used and proved to be stable after many tests.

• When dealing with "standard" log-likelihood minimization, the errors on the parameters are defined taking the value at the minimum + 0.5. This is at variance wrt the χ², where the errors are defined at the minimum + 1. Using a N-root of the likelihood translates into dinviding the log-likelihood by N, and the errors should be defined taking the value at the minimum + 0.5/N. Since in our case we use 4 roots with different N_i (number of hadrons for the 2 cells and 2 periods) , one possibility is to define the errors taking the minimum +0.5/<N_mean>, where <N_mean> is the arithmetic mean of the 4 N_i, namely the total number of hadrons divided by 4. (Note that if you are using 8 cells, you should divide by 8).
  By default MINUIT uses the errors definition for standard χ² minimization, so this behavior should be changed with a call:
  gMinuit->mnexcm( "SET ERR", arglist , 1, ierflg );
  where arglist[0] is 0.5/<N_mean>.

• By default MINUIT tries to achieve a fast convergence using the fewest possibile number of function calls, with the risk that the precision asked by the user is not obtained. This behavior can be avoided setting a flag that allows MINUIT to waste calls in order to be sure that the values found are precise. This is done with the command:
  gMinuit->Command("SET STRATEGY 2");
  

• Perform two calls of MIGRAD. In the first fit, the parameters are initialized with some numbers decided by the user (e.g. small number around 1-2 % for the asymmetries) The second fit uses as initial values the results obtained in the 1st fit. In case of large statistics and easy convergence, the results of the fits are the same. For some more critical cases, the second MIGRAD call is needed to have the best results.
  The default sets of MIGRAD have been changed as :
  • Step size set to 0.1 in the first MIGRAD fit and 0.01 in the second; This is done with the 4th parameter of the command gMinuit->mnparm(...), to be called before MIGRAD
  • Calls put at 5000, tolerance for EDM (Estimated Distance to Minimum) set to 0.5 (default is 0.1). This is done directly in the MIGRAD call, as:
    gMinuit->mnexcm( "MIGRAD", arglist, 2, ierflg );
    in which arglist[0] =5000 and arglist[1] = 0.5. If you set the parameters like this, you should see in the output something like:
    START MIGRAD MINIMIZATION. STRATEGY 2. CONVERGENCE WHEN EDM .LT. 5.00e-04
    In case the first fit does not converge, the fit is performed once again increasing the tolerance of a factor 4 (from 0.5 to 2, then from 2 to 8 and so on)
• Call to HESSE: this call is needed to have a better determination of the errors. The routine calculates the full matrix of second derivatives of the function wrt parameters, and inverts it.

Correlation coefficients

Minuit evaluates the full covariance matrix, that means also the correlation coefficients between the asymmetries are accessible. In order to access them , use the method gMinuit->mnemat(..) The diagonal terms are the squared errors and the off-diagonal terms (j,k) are the correlation coefficients times the errors of the corresponding asymmetries j and k.

Results with this UL version

The results on 2010 data for unidentified hadrons with the UL described above can be found here. -- FedericaSozzi - 2014-03-07