Re: Sorting a branch

Rene Brun (Rene.Brun@cern.ch)
Fri, 19 Jun 1998 09:02:59 +0200


Scott Sampson wrote:
>
> Hi,
>
> I have a Tree and am interested in a subset of event which pass some
> TCut. For those events which pass the TCut, I would like to take the
> values in the Branch X, and put them into some sort of array. (For
> simplicity let's assume X is a just a Float_t). Then I would like to sort
> this array. Is there an easy way to do all this? And do it fast? I've had
> trouble just trying to sort arrays.
>
> My goal is to be able to quickly calculate things like the following:
> For those events which satisfy TCut, what is the value X_cut, such that
> 95% of these events also satisfy X > X_cut?
>

You will find below an example of function to sort an array.
I translated this function to C++ starting from the Fortran
algorithm in Cernlib.
We will try to collect these kind of frequently used functions
like Sort, FindZeroOfFunction, Prob, Denlan in a TUtil class.

Rene Brun

void SortDown(Int_t n1, Float_t *a, Int_t *index, Bool_t down)
{
// sort the n1 elements of array a.
// In output the array index contains the indices of the sorted array.
// if down is false sort in increasing order (default is decreasing
order)
// This is a translation of the CERNLIB routine sortzv (M101)
// based on the quicksort algorithm

Int_t i,i1,n,i2,i3,i33,i222,iswap,n2;
Int_t i22 = 0;
Float_t ai;
n = n1;
if (n <= 0) return;
if (n == 1) {index[0] = 0; return;}
for (i=0;i<n;i++) index[i] = i+1;
for (i1=2;i1<=n;i1++) {
i3 = i1;
i33 = index[i3-1];
ai = a[i33-1];
while(1) {
i2 = i3/2;
if (i2 <= 0) break;
i22 = index[i2-1];
if (ai <= a[i22-1]) break;
index[i3-1] = i22;
i3 = i2;
}
index[i3-1] = i33;
}

while(1) {
i3 = index[n-1];
index[n-1] = index[0];
ai = a[i3-1];
n--;
if(n-1 < 0) {index[0] = i3; break;}
i1 = 1;
while(2) {
i2 = i1+i1;
if (i2 <= n) i22 = index[i2-1];
if (i2-n > 0) {index[i1-1] = i3; break;}
if (i2-n < 0) {
i222 = index[i2];
if (a[i22-1] - a[i222-1] < 0) {
i2++;
i22 = i222;
}
}
if (ai - a[i22-1] > 0) {index[i1-1] = i3; break;}
index[i1-1] = i22;
i1 = i2;
}
}
if (!down) return;
n2 = n1/2;
for (i=0;i<n1;i++) index[i]--;
for (i=0;i<n2;i++) {
iswap = index[i];
index[i] = index[n1-i-1];
index[n1-i-1] = iswap;
}
}