PMD/AliPMDClustering.cxx

   1 //-----------------------------------------------------//
   2 //                                                     //
   3 //  Source File : PMDClustering.cxx, Version 00        //
   4 //                                                     //
   5 //  Date   : September 26 2002                         //
   6 //                                                     //
   7 //  clustering code for alice pmd                      //
   8 //                                                     //
   9 //-----------------------------------------------------//
  10
  11 /*
  12    --------------------------------------------------------------------
  13    Code developed by S. C. Phatak, Institute of Physics,
  14    Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
  15    ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
  16    builds up superclusters and breaks them into clusters. The input is
  17    in array d[ndimx][ndimy] and cluster information is in array
  18    clusters[5][5000]. integer clno gives total number of clusters in the
  19    supermodule.
  20
  21    d, clno  and clusters are the only global ( public ) variables. Others
  22    are local ( private ) to the code.
  23
  24    At the moment, the data is read for whole detector ( all supermodules
  25    and pmd as well as cpv. This will have to be modify later )
  26
  27    LAST UPDATE  :  October 23, 2002
  28 -----------------------------------------------------------------------
  29 */
  30
  31
  32
  33 #include <TNtuple.h>
  34 #include <TObjArray.h>
  35 #include "AliPMDcluster.h"
  36 #include "AliPMDClustering.h"
  37 #include <stdio.h>
  38
  39 ClassImp(AliPMDClustering)
  40
  41 const double AliPMDClustering::pi=3.141593;
  42 const double AliPMDClustering::sqrth=0.8660254;  // sqrth = sqrt(3.)/2.
  43
  44
  45 AliPMDClustering::AliPMDClustering()
  46 {
  47   fMessage = 0;
  48   for(int i = 0; i < ndimx; i++)
  49     {
  50       for(int j = 0; j < ndimy; j++)
  51         {
  52           coord[0][i][j] = i+j/2.;
  53           coord[1][i][j] = sqrth*j;
  54         }
  55     }
  56 }
  57 AliPMDClustering::~AliPMDClustering()
  58 {
  59
  60 }
  61
  62 void AliPMDClustering::DoClust(int idet, int isup, double d1[72][72], TObjArray *pmdcont)
  63 {
  64
  65   AliPMDcluster *pmdcl = 0;
  66
  67   int i, i1, i2, j, nmx1, incr;
  68   double  cutoff, ave;
  69   Float_t clusdata[5];
  70
  71   const float twobysqrt3 = 1.1547; // 2./sqrt(3.)
  72
  73   //  if (fMessage == 1)
  74     {
  75       cout << " supermodule no. " << idet << " " << isup << endl;
  76     }
  77
  78   for (i = 0; i < ndimx; i++)
  79     {
  80       for (j = 0; j < ndimy; j++)
  81         {
  82           d[i][j] = d1[i][j];
  83         }
  84     }
  85   order(idet); // order the data
  86   cutoff=400.; // cutoff used to discard cells having ener. dep.
  87   ave=0.;
  88   nmx1=-1;
  89
  90   for(j=0;j<nmx; j++)
  91     {
  92       i1 = iord[0][j];
  93       i2 = iord[1][j];
  94       if (d[i1][i2] > 0.) {ave=ave+d[i1][i2];}
  95       if (d[i1][i2] >= cutoff ) nmx1 = nmx1 + 1;
  96     }
  97   // nmx1 --- number of cells having ener dep >= cutoff
  98   if (fMessage == 1)
  99     {
 100       cout << " nmx1 " << nmx1 << endl;
 101     }
 102   ave=ave/nmx1;
 103   if (fMessage == 1)
 104     {
 105       cout <<"nmx " << nmx << " nmx1 " << nmx1<< " ave "<<ave<<
 106         " cutoff " << cutoff << endl;
 107     }
 108
 109   incr = crclust(ave, cutoff, nmx1, idet);
 110
 111   refclust(incr, i, idet);
 112   if (fMessage == 1)
 113     {
 114       if(idet == 0)cout <<" supermodule ( pmd ) = "<< isup <<" done "
 115                         <<endl;
 116       if(idet == 1)cout <<" supermodule ( cpv ) = "<< isup <<" done "
 117                         <<endl;
 118       cout << "clno " << clno << endl;
 119     }
 120
 121
 122   for(i1=0; i1<clno; i1++)
 123     {
 124       float clu_xc    = (float) clusters[0][i1];
 125       float clu_yc    = (float) clusters[1][i1];
 126       float clu_adc   = (float) clusters[2][i1];
 127       float clu_cells = (float) clusters[3][i1];
 128       float clu_rad   = (float) clusters[4][i1];
 129
 130       float clu_y0 = twobysqrt3*clu_yc;
 131       float clu_x0 = clu_xc - clu_y0/2.;
 132
 133       clusdata[0] = clu_cells;
 134       clusdata[1] = clu_x0;
 135       clusdata[2] = clu_y0;
 136       clusdata[3] = clu_adc;
 137       clusdata[4] = clu_rad;
 138
 139       pmdcl = new AliPMDcluster(clusdata);
 140       pmdcont->Add(pmdcl);
 141     }
 142
 143   delete pmdcl;
 144
 145 }
 146
 147 void AliPMDClustering::order(int /*idet*/)
 148 {
 149   // using simple sort
 150   double dd[nmx], adum;// matrix d converted into
 151   // one dimensional array dd. adum a place holder for double
 152   int i, j, i1, i2, iord1[nmx], itst, idum; // information of
 153   // ordering is stored in iord1, original array not ordered
 154   //
 155   // define arrays dd and iord1
 156   for(i1=0; i1 < ndimx; i1++){
 157     for(i2=0; i2 < ndimy; i2++){
 158       i=i1+i2*ndimx;
 159       iord1[i]=i; dd[i]=d[i1][i2];
 160     }
 161   }
 162   // sort and store sorting information in iord1
 163   for(j=1; j < nmx; j++){
 164     itst=0; adum=dd[j]; idum=iord1[j];
 165     for(i1=0; i1 < j ; i1++){
 166       if(adum > dd[i1] && itst == 0){
 167         itst=1;
 168         for(i2=j-1; i2 >= i1 ; i2=i2--){
 169           dd[i2+1]=dd[i2];
 170           iord1[i2+1]=iord1[i2];
 171         }
 172         dd[i1]=adum; iord1[i1]=idum;
 173       }
 174     }
 175   }
 176   // store the sorted information in iord for later use
 177   for(i=0; i<nmx; i++){
 178     j=iord1[i]; i2=j/ndimx;
 179     i1=j-i2*ndimx;
 180     iord[0][i]=i1;
 181     iord[1][i]=i2;
 182   }
 183 }
 184
 185 void AliPMDClustering::refclust(int incr, int /*supmod*/, int /*idet*/)
 186 {
 187   int i, j, k, i1, i2, id, icl, ncl[4500], iord[4500], itest;
 188   int ihld;
 189   int ig, nsupcl, lev1[20], lev2[20];
 190   double x[4500], y[4500], z[4500], x1, y1, z1, x2, y2, z2, dist;
 191   double xc[4500], yc[4500], zc[4500], cells[4500], sum, rc[4500], rr;
 192   // clno counts the final clusters
 193   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 194   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 195   // xc, yc store (x,y) coordinates of the cluster center
 196   // zc stores the energy deposited in a cluster
 197   // rc is cluster radius
 198   // finally the cluster information is put in 2-dimensional array clusters
 199   //  ofstream ofl1("checking.5",ios::app);
 200   clno=-1;
 201   nsupcl=-1;
 202   for(i=0; i<4500; i++){ncl[i]=-1;}
 203   for(i=0; i<incr; i++){
 204     if(infcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
 205     ncl[nsupcl]=ncl[nsupcl]+1;
 206   }
 207   if (fMessage == 1)
 208     {
 209       cout << " # of cells " <<incr+1 << " # of superclusters " << nsupcl+1
 210            << endl;
 211     }
 212   id=-1;
 213   icl=-1;
 214   for(i=0; i<nsupcl; i++){
 215     if(ncl[i] == 0){
 216       id=id+1;
 217       icl=icl+1;
 218       // one  cell super-clusters --> single cluster
 219       // cluster center at the centyer of the cell
 220       // cluster radius = half cell dimension
 221       clno=clno+1;
 222       i1=infcl[1][id];
 223       i2=infcl[2][id];
 224       clusters[0][clno]=coord[0][i1][i2];
 225       clusters[1][clno]=coord[1][i1][i2];
 226       clusters[2][clno]=d[i1][i2];
 227       clusters[3][clno]=1.;
 228       clusters[4][clno]=0.5;
 229       //ofl1 << icl << " " << coord[0][i1][i2] << " " << coord[1][i1][i2] <<
 230       //" " << d[i1][i2] << " " << clusters[3][clno] <<endl;
 231     }else if(ncl[i] == 1){
 232       // two cell super-cluster --> single cluster
 233       // cluster center is at ener. dep.-weighted mean of two cells
 234       // cluster radius == half cell dimension
 235       id=id+1;
 236       icl=icl+1;
 237       clno=clno+1;
 238       i1=infcl[1][id];
 239       i2=infcl[2][id];
 240       x1=coord[0][i1][i2];
 241       y1=coord[1][i1][i2];
 242       z1=d[i1][i2];
 243       id=id+1;
 244       i1=infcl[1][id];
 245       i2=infcl[2][id];
 246       x2=coord[0][i1][i2];
 247       y2=coord[1][i1][i2];
 248       z2=d[i1][i2];
 249       clusters[0][clno]=(x1*z1+x2*z2)/(z1+z2);
 250       clusters[1][clno]=(y1*z1+y2*z2)/(z1+z2);
 251       clusters[2][clno]=z1+z2;
 252       clusters[3][clno]=2.;
 253       clusters[4][clno]=0.5;
 254       //ofl1 << icl << " " << clusters[0][clno] << " " << clusters[1][clno]
 255       //   << " " << clusters[2][clno] << " " <<clusters[3][clno] <<endl;
 256     }else{
 257       id=id+1;
 258       iord[0]=0;
 259       // super-cluster of more than two cells - broken up into smaller
 260       // clusters gaussian centers computed. (peaks separated by > 1 cell)
 261       // Begin from cell having largest energy deposited This is first
 262       // cluster center
 263       i1=infcl[1][id];
 264       i2=infcl[2][id];
 265       x[0]=coord[0][i1][i2];
 266       y[0]=coord[1][i1][i2];
 267       z[0]=d[i1][i2];
 268       iord[0]=0;
 269       for(j=1;j<=ncl[i];j++){
 270         id=id+1;
 271         i1=infcl[1][id];
 272         i2=infcl[2][id];
 273         iord[j]=j;
 274         x[j]=coord[0][i1][i2];
 275         y[j]=coord[1][i1][i2];
 276         z[j]=d[i1][i2];
 277       }
 278       // arranging cells within supercluster in decreasing order
 279       for(j=1;j<=ncl[i];j++){
 280         itest=0; ihld=iord[j];
 281         for(i1=0;i1<j;i1++){
 282           if(itest == 0 && z[iord[i1]] < z[ihld]){
 283             itest=1;
 284             for(i2=j-1;i2>=i1;i2--){
 285               iord[i2+1]=iord[i2];
 286             }
 287             iord[i1]=ihld;
 288           }
 289         }
 290       }
 291       // compute the number of Gaussians and their centers ( first
 292       // guess )
 293       // centers must be separated by cells having smaller ener. dep.
 294       // neighbouring centers should be either strong or well-separated
 295       ig=0;
 296       xc[ig]=x[iord[0]];
 297       yc[ig]=y[iord[0]];
 298       zc[ig]=z[iord[0]];
 299       for(j=1;j<=ncl[i];j++){
 300         itest=-1;
 301         x1=x[iord[j]];
 302         y1=y[iord[j]];
 303         for(k=0;k<=ig;k++){
 304           x2=xc[k]; y2=yc[k];
 305           rr=Dist(x1,y1,x2,y2);
 306           if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
 307             itest=itest+1;
 308           if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
 309             itest=itest+1;
 310           if( rr >= 2.1)itest=itest+1;
 311         }
 312         if(itest == ig){
 313           ig=ig+1;
 314           xc[ig]=x1;
 315           yc[ig]=y1;
 316           zc[ig]=z[iord[j]];
 317         }
 318       }
 319       // for(j=0; j<=ig; j++){
 320       //ofl1 << icl+j+1 << " " << xc[j] << " " <<yc[j] <<" "<<zc[j]<<endl;
 321       //}
 322       // gaussfit to adjust cluster parameters to minimize
 323       gaussfit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
 324       icl=icl+ig+1;
 325       // compute the number of cells belonging to each cluster.
 326       // cell is shared between several clusters ( if they are equidistant
 327       // from it ) in the ratio of cluster energy deposition
 328       for(j=0; j<=ig; j++){
 329         cells[j]=0.;
 330       }
 331       if(ig > 0){
 332         for(j=0; j<=ncl[i]; j++){
 333           lev1[0]=0;
 334           lev2[0]=0;
 335           for(k=0; k<=ig; k++){
 336             dist=Dist(x[j], y[j], xc[k], yc[k]);
 337             if(dist < sqrt(3.) ){
 338               lev1[0]++;
 339               i1=lev1[0];
 340               lev1[i1]=k;
 341             }else{
 342               if(dist < 2.1){
 343                 lev2[0]++;
 344                 i1=lev2[0];
 345                 lev2[i1]=k;
 346               }
 347             }
 348           }
 349           if(lev1[0] != 0){
 350             if(lev1[0] == 1){cells[lev1[1]]=cells[lev1[1]]+1.;}
 351             else{
 352               sum=0.;
 353               for(k=1; k<=lev1[0]; k++){
 354                 sum=sum+zc[lev1[k]];
 355               }
 356               for(k=1; k<=lev1[0]; k++){
 357                 cells[lev1[k]]=cells[lev1[k]]+zc[lev1[k]]/sum;
 358               }
 359             }
 360           }else{
 361             if(lev2[0] == 0){cells[lev2[1]]=cells[lev2[1]]+1.;}
 362             else{
 363               sum=0.;
 364               for(k=1; k<=lev2[0]; k++){
 365                 sum=sum+zc[lev2[k]];
 366               }
 367               for(k=1; k<=lev2[0]; k++){
 368                 cells[lev2[k]]=cells[lev2[k]]+zc[lev2[k]]/sum;
 369               }
 370             }
 371           }
 372         }
 373       }
 374       for(j=0; j<=ig; j++){
 375         clno=clno+1;
 376         clusters[0][clno]=xc[j];
 377         clusters[1][clno]=yc[j];
 378         clusters[2][clno]=zc[j];
 379         clusters[4][clno]=rc[j];
 380         if(ig == 0){
 381           clusters[3][clno]=ncl[i];
 382         }else{
 383           clusters[3][clno]=cells[j];
 384         }
 385       }
 386     }
 387   }
 388
 389 }
 390
 391
 392 void AliPMDClustering::gaussfit(int ncell, int nclust, double &x, double &y ,double &z, double &xc, double &yc, double &zc, double &rc)
 393 {
 394   int i, j, i1, i2, jmax, novar, idd, jj;
 395   double xx[4500], yy[4500], zz[4500], xxc[4500], yyc[4500];
 396   double a[4500], b[4500], c[4500], d[4500], ha[4500], hb[4500];
 397   double hc[4500], hd[4500], zzc[4500], rrc[4500];
 398   int neib[4500][50];
 399   double sum, dx, dy, str, str1, aint, sum1, rr, dum;
 400   double x1, x2, y1, y2;
 401   str=0.;
 402   str1=0.;
 403   rr=0.3;
 404   novar=0;
 405
 406   j = 0;  // Just put not to see the compiler warning, BKN
 407
 408
 409   for(i=0; i<=ncell; i++){
 410     xx[i]=*(&x+i);
 411     yy[i]=*(&y+i);
 412     zz[i]=*(&z+i);
 413     str=str+zz[i];
 414   }
 415   for(i=0; i<=nclust; i++){
 416     xxc[i]=*(&xc+i);
 417     yyc[i]=*(&yc+i);
 418     zzc[i]=*(&zc+i);
 419     str1=str1+zzc[i];
 420     rrc[i]=0.5;
 421
 422   }
 423   for(i=0; i<=nclust; i++){
 424     zzc[i]=str/str1*zzc[i];
 425     ha[i]=xxc[i];
 426     hb[i]=yyc[i];
 427     hc[i]=zzc[i];
 428     hd[i]=rrc[i];
 429     x1=xxc[i];
 430     y1=yyc[i];
 431   }
 432   for(i=0; i<=ncell; i++){
 433     idd=0;
 434     x1=xx[i];
 435     y1=yy[i];
 436     for(j=0; j<=nclust; j++){
 437       x2=xxc[j];
 438       y2=yyc[j];
 439       if(Dist(x1,y1,x2,y2) <= 3.){ idd=idd+1; neib[i][idd]=j; }
 440     }
 441
 442     neib[i][0]=idd;
 443   }
 444   sum=0.;
 445   for(i1=0; i1<=ncell; i1++){
 446     aint=0.;
 447     idd=neib[i1][0];
 448     for(i2=1; i2<=idd; i2++){
 449       jj=neib[i1][i2];
 450       dx=xx[i1]-xxc[jj];
 451       dy=yy[i1]-yyc[jj];
 452       dum=rrc[j]*rrc[jj]+rr*rr;
 453       aint=aint+exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
 454     }
 455     sum=sum+(aint-zz[i1])*(aint-zz[i1])/str;
 456   }
 457   jmax=nclust*1000;
 458   if(nclust > 20)jmax=20000;
 459   for(j=0; j<jmax; j++){
 460     str1=0.;
 461     for(i=0; i<=nclust; i++){
 462       a[i]=xxc[i]+0.6*(ranmar()-0.5);
 463       b[i]=yyc[i]+0.6*(ranmar()-0.5);
 464       c[i]=zzc[i]*(1.+(ranmar()-0.5)*0.2);
 465       str1=str1+zzc[i];
 466       d[i]=rrc[i]*(1.+(ranmar()-0.5)*0.1);
 467       if(d[i] < 0.25)d[i]=0.25;
 468     }
 469     for(i=0; i<=nclust; i++){ c[i]=c[i]*str/str1; }
 470     sum1=0.;
 471     for(i1=0; i1<=ncell; i1++){
 472       aint=0.;
 473       idd=neib[i1][0];
 474       for(i2=1; i2<=idd; i2++){
 475         jj=neib[i1][i2];
 476         dx=xx[i1]-a[jj];
 477         dy=yy[i1]-b[jj];
 478         dum=d[jj]*d[jj]+rr*rr;
 479         aint=aint+exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
 480       }
 481       sum1=sum1+(aint-zz[i1])*(aint-zz[i1])/str;
 482     }
 483
 484     if(sum1 < sum){
 485       for(i2=0; i2<=nclust; i2++){
 486         xxc[i2]=a[i2];
 487         yyc[i2]=b[i2];
 488         zzc[i2]=c[i2];
 489         rrc[i2]=d[i2];
 490         sum=sum1;
 491
 492       }
 493     }
 494   }
 495   for(j=0; j<=nclust; j++){
 496     *(&xc+j)=xxc[j];
 497     *(&yc+j)=yyc[j];
 498     *(&zc+j)=zzc[j];
 499     *(&rc+j)=rrc[j];
 500   }
 501 }
 502
 503
 504 double AliPMDClustering::Dist(double x1, double y1, double x2, double y2)
 505 {
 506   return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 507 }
 508
 509
 510 int AliPMDClustering::crclust(double /*ave*/, double cutoff, int nmx1, int /*idet*/)
 511 {
 512   int i,j,k,id1,id2,icl, numcell, clust[2][5000];
 513   int jd1,jd2, icell, cellcount;
 514   static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 515   // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
 516   // cell. There are six neighbours.
 517   // cellcount --- total number of cells having nonzero ener dep
 518   // numcell --- number of cells in a given supercluster
 519   //ofstream ofl0("cells_loc",ios::out);
 520   // initialize infocl[2][ndimx][ndimy]
 521   for (j=0; j < 72; j++){
 522     for(k=0; k < 72; k++){
 523       infocl[0][j][k] = 0;
 524       infocl[1][j][k] = 0;
 525     }
 526   }
 527   for(i=0; i < nmx; i++){
 528     infcl[0][i] = -1;
 529     id1=iord[0][i];
 530     id2=iord[1][i];
 531     if(d[id1][id2] <= cutoff){infocl[0][id1][id2]=-1;}
 532   }
 533   // ---------------------------------------------------------------
 534   // crude clustering begins. Start with cell having largest adc
 535   // count and loop over the cells in descending order of adc count
 536   // ---------------------------------------------------------------
 537   icl=-1;
 538   cellcount=-1;
 539   for(icell=0; icell <= nmx1; icell++){
 540     id1=iord[0][icell];
 541     id2=iord[1][icell];
 542     if(infocl[0][id1][id2] == 0 ){
 543       // ---------------------------------------------------------------
 544       // icl -- cluster #, numcell -- # of cells in it, clust -- stores
 545       // coordinates of the cells in a cluster, infocl[0][i1][i2] is 1 for
 546       // primary and 2 for secondary cells,
 547       // infocl[1][i1][i2] stores cluster #
 548       // ---------------------------------------------------------------
 549       icl=icl+1;
 550       numcell=0;
 551       cellcount=cellcount+1;
 552       infocl[0][id1][id2]=1;
 553       infocl[1][id1][id2]=icl;
 554       infcl[0][cellcount]=icl;
 555       infcl[1][cellcount]=id1;
 556       infcl[2][cellcount]=id2;
 557       clust[0][numcell]=id1;
 558       clust[1][numcell]=id2;
 559       for(i=1; i<5000; i++)clust[0][i]=0;
 560       // ---------------------------------------------------------------
 561       // check for adc count in neib. cells. If ne 0 put it in this clust
 562       // ---------------------------------------------------------------
 563       for(i=0; i<6; i++){
 564         jd1=id1+neibx[i];
 565         jd2=id2+neiby[i];
 566         if( (jd1 >= 0 && jd1 < 72) && (jd2 >= 0 && jd2 < 72) &&
 567             infocl[0][jd1][jd2] == 0){
 568           numcell=numcell+1;
 569           infocl[0][jd1][jd2]=2;
 570           infocl[1][jd1][jd2]=icl;
 571           clust[0][numcell]=jd1;
 572           clust[1][numcell]=jd2;
 573           cellcount=cellcount+1;
 574           infcl[0][cellcount]=icl;
 575           infcl[1][cellcount]=jd1;
 576           infcl[2][cellcount]=jd2;
 577         }
 578       }
 579       // ---------------------------------------------------------------
 580       // check adc count for neighbour's neighbours recursively and
 581       // if nonzero, add these to the cluster.
 582       // ---------------------------------------------------------------
 583       for(i=1;i < 5000;i++){
 584         if(clust[0][i] != 0){
 585           id1=clust[0][i];
 586           id2=clust[1][i];
 587           for(j=0; j<6 ; j++){
 588             jd1=id1+neibx[j];
 589             jd2=id2+neiby[j];
 590             if( (jd1 >= 0 && jd1 < 72) && (jd2 >= 0 && jd2 < 72) &&
 591                 infocl[0][jd1][jd2] == 0 ){
 592               infocl[0][jd1][jd2]=2;
 593               infocl[1][jd1][jd2]=icl;
 594               numcell=numcell+1;
 595               clust[0][numcell]=jd1;
 596               clust[1][numcell]=jd2;
 597               cellcount=cellcount+1;
 598               infcl[0][cellcount]=icl;
 599               infcl[1][cellcount]=jd1;
 600               infcl[2][cellcount]=jd2;
 601             }
 602           }
 603         }
 604       }
 605     }
 606   }
 607   //  for(icell=0; icell<=cellcount; icell++){
 608   //    ofl0 << infcl[0][icell] << " " << infcl[1][icell] << " " <<
 609   //      infcl[2][icell] << endl;
 610   //  }
 611   return cellcount;
 612 }
 613
 614 double AliPMDClustering::ranmar()
 615 {
 616   /*                                   C==========================C*/
 617   /*===================================C==========================*/
 618   /*  Universal random number generator proposed by Marsaglia and Zaman
 619       in report FSU-SCRI-87-50 */
 620
 621   //  clock_t start;
 622   int ii, jj;
 623   static int i=96, j=32, itest=0, i1, i2, i3, i4, i5;
 624   static double u[97], c, cd, cm, s, t;
 625   static double uni;
 626   int count1,count2,idum;
 627   /*    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  */
 628   if (itest == 0) {
 629     //*******************************************************
 630     // following three lines if the seed to be provided by computer
 631     // start = time(NULL);
 632     // ii=start;
 633     // jj=start;
 634     //*******************************************************
 635     //following two lines for fixed seed ( during testing only. Else
 636     //use preceeing three lines
 637     ii=8263;
 638     jj=5726;
 639     if(ii > 31328 ) ii = ii - ( ii / 31328 ) * 31328;
 640     if(jj > 30081 ) jj = jj - ( jj / 30081 ) * 30081;
 641     itest=itest+1;
 642     if((( ii > 0 ) &&  ( ii <= 31328 )) && (( jj > 0 ) &&
 643                                             ( jj <= 30081 ))){
 644       i1=ii/177+2; i2=ii-(i1-2)*177+2; i3=jj/169+1; i4=jj-(i3-1)*169;
 645       i4 = jj - (i3-1)*169;
 646       count1=0;
 647       while ( count1 < 97 ){
 648         s=0.;
 649         t=0.5;
 650         count2=0;
 651         while( count2 < 24 ){
 652           idum=i1*i2/179;
 653           idum=( i1*i2 - (i1*i2/179)*179 ) * i3;
 654           i5=idum-(idum/179)*179;
 655           i1=i2; i2=i3; i3=i5; idum=53*i4+1; i4=idum-(idum/169)*169;
 656           if( i4*i5-((i4*i5)/64)*64 >= 32 ) s=s+t;
 657           t=0.5*t;
 658           count2=count2+1;
 659         }
 660         u[count1] = s;
 661         count1 = count1 +1;
 662       }
 663       c = 362436./16777216.;  cd = 7654321./16777216.;
 664       cm = 16777213./16777216.;
 665     }
 666     else{
 667       cout << " wrong initialization " << endl;
 668     }
 669   }
 670   else{
 671     uni = u[i] - u[j]; if( uni < 0.) uni = uni + 1; u[i] = uni;
 672     i = i -1;
 673     if( i < 0 ) i = 96; j = j - 1; if ( j < 0 ) j = 96; c = c - cd;
 674     if( c < 0. ) c = c+cm; uni = uni-c ; if( uni < 0. )uni = uni+1.;
 675     //    return uni;
 676   }
 677   return uni;
 678
 679 }
 680
 681 void AliPMDClustering::SetMessage(Int_t imessage)
 682 {
 683   fMessage = imessage;
 684 }