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 AliPMDClustering::AliPMDClustering()
  42 {
  43   fMessage = 0;
  44   for(int i = 0; i < ndimx; i++)
  45     {
  46       for(int j = 0; j < ndimy; j++)
  47         {
  48           coord[0][i][j] = i+j/2.;
  49           coord[1][i][j] = sqrth*j;
  50         }
  51     }
  52 }
  53 AliPMDClustering::~AliPMDClustering()
  54 {
  55
  56 }
  57
  58 //void AliPMDClustering::DoClust(int idet, int isup, double d1[72][72], AliPMDContainer *pmdc)
  59 void AliPMDClustering::DoClust(int idet, int isup, double d1[72][72], TObjArray *pmdcont)
  60 {
  61
  62   AliPMDcluster *pmdcl = 0;
  63
  64   int i, i1, i2, j, nmx1, incr;
  65   double  cutoff, ave;
  66   Float_t clusdata[5];
  67
  68   const float twobysqrt3 = 1.1547; // 2./sqrt(3.)
  69
  70   //  if (fMessage == 1)
  71     {
  72       cout << " supermodule no. " << idet << " " << isup << endl;
  73     }
  74
  75   for (i = 0; i < ndimx; i++)
  76     {
  77       for (j = 0; j < ndimy; j++)
  78         {
  79           d[i][j] = d1[i][j];
  80         }
  81     }
  82   order(idet); // order the data
  83   cutoff=400.; // cutoff used to discard cells having ener. dep.
  84   ave=0.;
  85   nmx1=-1;
  86
  87   for(j=0;j<nmx; j++)
  88     {
  89       i1 = iord[0][j];
  90       i2 = iord[1][j];
  91       if (d[i1][i2] > 0.) {ave=ave+d[i1][i2];}
  92       if (d[i1][i2] >= cutoff ) nmx1 = nmx1 + 1;
  93     }
  94   // nmx1 --- number of cells having ener dep >= cutoff
  95   if (fMessage == 1)
  96     {
  97       cout << " nmx1 " << nmx1 << endl;
  98     }
  99   ave=ave/nmx1;
 100   if (fMessage == 1)
 101     {
 102       cout <<"nmx " << nmx << " nmx1 " << nmx1<< " ave "<<ave<<
 103         " cutoff " << cutoff << endl;
 104     }
 105
 106   incr = crclust(ave, cutoff, nmx1, idet);
 107
 108   refclust(incr, i, idet);
 109   if (fMessage == 1)
 110     {
 111       if(idet == 0)cout <<" supermodule ( pmd ) = "<< isup <<" done "
 112                         <<endl;
 113       if(idet == 1)cout <<" supermodule ( cpv ) = "<< isup <<" done "
 114                         <<endl;
 115       cout << "clno " << clno << endl;
 116     }
 117
 118
 119   for(i1=0; i1<clno; i1++)
 120     {
 121       float clu_xc    = (float) clusters[0][i1];
 122       float clu_yc    = (float) clusters[1][i1];
 123       float clu_adc   = (float) clusters[2][i1];
 124       float clu_cells = (float) clusters[3][i1];
 125       float clu_rad   = (float) clusters[4][i1];
 126
 127       float clu_y0 = twobysqrt3*clu_yc;
 128       float clu_x0 = clu_xc - clu_y0/2.;
 129
 130       clusdata[0] = clu_cells;
 131       clusdata[1] = clu_x0;
 132       clusdata[2] = clu_y0;
 133       clusdata[3] = clu_adc;
 134       clusdata[4] = clu_rad;
 135
 136       pmdcl = new AliPMDcluster(clusdata);
 137       pmdcont->Add(pmdcl);
 138     }
 139
 140   delete pmdcl;
 141
 142 }
 143
 144 void AliPMDClustering::order(int idet)
 145 {
 146   // using simple sort
 147   double dd[nmx], adum;// matrix d converted into
 148   // one dimensional array dd. adum a place holder for double
 149   int i, j, i1, i2, iord1[nmx], itst, idum; // information of
 150   // ordering is stored in iord1, original array not ordered
 151   //
 152   // define arrays dd and iord1
 153   for(i1=0; i1 < ndimx; i1++){
 154     for(i2=0; i2 < ndimy; i2++){
 155       i=i1+i2*ndimx;
 156       iord1[i]=i; dd[i]=d[i1][i2];
 157     }
 158   }
 159   // sort and store sorting information in iord1
 160   for(j=1; j < nmx; j++){
 161     itst=0; adum=dd[j]; idum=iord1[j];
 162     for(i1=0; i1 < j ; i1++){
 163       if(adum > dd[i1] && itst == 0){
 164         itst=1;
 165         for(i2=j-1; i2 >= i1 ; i2=i2--){
 166           dd[i2+1]=dd[i2];
 167           iord1[i2+1]=iord1[i2];
 168         }
 169         dd[i1]=adum; iord1[i1]=idum;
 170       }
 171     }
 172   }
 173   // store the sorted information in iord for later use
 174   for(i=0; i<nmx; i++){
 175     j=iord1[i]; i2=j/ndimx;
 176     i1=j-i2*ndimx;
 177     iord[0][i]=i1;
 178     iord[1][i]=i2;
 179   }
 180 }
 181
 182 void AliPMDClustering::refclust(int incr, int supmod, int idet)
 183 {
 184   int i, j, k, i1, i2, id, icl, ncl[4500], iord[4500], itest;
 185   int ihld;
 186   int ig, nsupcl, lev1[20], lev2[20];
 187   double x[4500], y[4500], z[4500], x1, y1, z1, x2, y2, z2, dist;
 188   double xc[4500], yc[4500], zc[4500], cells[4500], sum, rc[4500], rr;
 189   // clno counts the final clusters
 190   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 191   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 192   // xc, yc store (x,y) coordinates of the cluster center
 193   // zc stores the energy deposited in a cluster
 194   // rc is cluster radius
 195   // finally the cluster information is put in 2-dimensional array clusters
 196   //  ofstream ofl1("checking.5",ios::app);
 197   clno=-1;
 198   nsupcl=-1;
 199   for(i=0; i<4500; i++){ncl[i]=-1;}
 200   for(i=0; i<incr; i++){
 201     if(infcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
 202     ncl[nsupcl]=ncl[nsupcl]+1;
 203   }
 204   if (fMessage == 1)
 205     {
 206       cout << " # of cells " <<incr+1 << " # of superclusters " << nsupcl+1
 207            << endl;
 208     }
 209   id=-1;
 210   icl=-1;
 211   for(i=0; i<nsupcl; i++){
 212     if(ncl[i] == 0){
 213       id=id+1;
 214       icl=icl+1;
 215       // one  cell super-clusters --> single cluster
 216       // cluster center at the centyer of the cell
 217       // cluster radius = half cell dimension
 218       clno=clno+1;
 219       i1=infcl[1][id];
 220       i2=infcl[2][id];
 221       clusters[0][clno]=coord[0][i1][i2];
 222       clusters[1][clno]=coord[1][i1][i2];
 223       clusters[2][clno]=d[i1][i2];
 224       clusters[3][clno]=1.;
 225       clusters[4][clno]=0.5;
 226       //ofl1 << icl << " " << coord[0][i1][i2] << " " << coord[1][i1][i2] <<
 227       //" " << d[i1][i2] << " " << clusters[3][clno] <<endl;
 228     }else if(ncl[i] == 1){
 229       // two cell super-cluster --> single cluster
 230       // cluster center is at ener. dep.-weighted mean of two cells
 231       // cluster radius == half cell dimension
 232       id=id+1;
 233       icl=icl+1;
 234       clno=clno+1;
 235       i1=infcl[1][id];
 236       i2=infcl[2][id];
 237       x1=coord[0][i1][i2];
 238       y1=coord[1][i1][i2];
 239       z1=d[i1][i2];
 240       id=id+1;
 241       i1=infcl[1][id];
 242       i2=infcl[2][id];
 243       x2=coord[0][i1][i2];
 244       y2=coord[1][i1][i2];
 245       z2=d[i1][i2];
 246       clusters[0][clno]=(x1*z1+x2*z2)/(z1+z2);
 247       clusters[1][clno]=(y1*z1+y2*z2)/(z1+z2);
 248       clusters[2][clno]=z1+z2;
 249       clusters[3][clno]=2.;
 250       clusters[4][clno]=0.5;
 251       //ofl1 << icl << " " << clusters[0][clno] << " " << clusters[1][clno]
 252       //   << " " << clusters[2][clno] << " " <<clusters[3][clno] <<endl;
 253     }else{
 254       id=id+1;
 255       iord[0]=0;
 256       // super-cluster of more than two cells - broken up into smaller
 257       // clusters gaussian centers computed. (peaks separated by > 1 cell)
 258       // Begin from cell having largest energy deposited This is first
 259       // cluster center
 260       i1=infcl[1][id];
 261       i2=infcl[2][id];
 262       x[0]=coord[0][i1][i2];
 263       y[0]=coord[1][i1][i2];
 264       z[0]=d[i1][i2];
 265       iord[0]=0;
 266       for(j=1;j<=ncl[i];j++){
 267         id=id+1;
 268         i1=infcl[1][id];
 269         i2=infcl[2][id];
 270         iord[j]=j;
 271         x[j]=coord[0][i1][i2];
 272         y[j]=coord[1][i1][i2];
 273         z[j]=d[i1][i2];
 274       }
 275       // arranging cells within supercluster in decreasing order
 276       for(j=1;j<=ncl[i];j++){
 277         itest=0; ihld=iord[j];
 278         for(i1=0;i1<j;i1++){
 279           if(itest == 0 && z[iord[i1]] < z[ihld]){
 280             itest=1;
 281             for(i2=j-1;i2>=i1;i2--){
 282               iord[i2+1]=iord[i2];
 283             }
 284             iord[i1]=ihld;
 285           }
 286         }
 287       }
 288       // compute the number of Gaussians and their centers ( first
 289       // guess )
 290       // centers must be separated by cells having smaller ener. dep.
 291       // neighbouring centers should be either strong or well-separated
 292       ig=0;
 293       xc[ig]=x[iord[0]];
 294       yc[ig]=y[iord[0]];
 295       zc[ig]=z[iord[0]];
 296       for(j=1;j<=ncl[i];j++){
 297         itest=-1;
 298         x1=x[iord[j]];
 299         y1=y[iord[j]];
 300         for(k=0;k<=ig;k++){
 301           x2=xc[k]; y2=yc[k];
 302           rr=Dist(x1,y1,x2,y2);
 303           if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
 304             itest=itest+1;
 305           if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
 306             itest=itest+1;
 307           if( rr >= 2.1)itest=itest+1;
 308         }
 309         if(itest == ig){
 310           ig=ig+1;
 311           xc[ig]=x1;
 312           yc[ig]=y1;
 313           zc[ig]=z[iord[j]];
 314         }
 315       }
 316       // for(j=0; j<=ig; j++){
 317       //ofl1 << icl+j+1 << " " << xc[j] << " " <<yc[j] <<" "<<zc[j]<<endl;
 318       //}
 319       // gaussfit to adjust cluster parameters to minimize
 320       gaussfit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
 321       icl=icl+ig+1;
 322       // compute the number of cells belonging to each cluster.
 323       // cell is shared between several clusters ( if they are equidistant
 324       // from it ) in the ratio of cluster energy deposition
 325       for(j=0; j<=ig; j++){
 326         cells[j]=0.;
 327       }
 328       if(ig > 0){
 329         for(j=0; j<=ncl[i]; j++){
 330           lev1[0]=0;
 331           lev2[0]=0;
 332           for(k=0; k<=ig; k++){
 333             dist=Dist(x[j], y[j], xc[k], yc[k]);
 334             if(dist < sqrt(3.) ){
 335               lev1[0]++;
 336               i1=lev1[0];
 337               lev1[i1]=k;
 338             }else{
 339               if(dist < 2.1){
 340                 lev2[0]++;
 341                 i1=lev2[0];
 342                 lev2[i1]=k;
 343               }
 344             }
 345           }
 346           if(lev1[0] != 0){
 347             if(lev1[0] == 1){cells[lev1[1]]=cells[lev1[1]]+1.;}
 348             else{
 349               sum=0.;
 350               for(k=1; k<=lev1[0]; k++){
 351                 sum=sum+zc[lev1[k]];
 352               }
 353               for(k=1; k<=lev1[0]; k++){
 354                 cells[lev1[k]]=cells[lev1[k]]+zc[lev1[k]]/sum;
 355               }
 356             }
 357           }else{
 358             if(lev2[0] == 0){cells[lev2[1]]=cells[lev2[1]]+1.;}
 359             else{
 360               sum=0.;
 361               for(k=1; k<=lev2[0]; k++){
 362                 sum=sum+zc[lev2[k]];
 363               }
 364               for(k=1; k<=lev2[0]; k++){
 365                 cells[lev2[k]]=cells[lev2[k]]+zc[lev2[k]]/sum;
 366               }
 367             }
 368           }
 369         }
 370       }
 371       for(j=0; j<=ig; j++){
 372         clno=clno+1;
 373         clusters[0][clno]=xc[j];
 374         clusters[1][clno]=yc[j];
 375         clusters[2][clno]=zc[j];
 376         clusters[4][clno]=rc[j];
 377         if(ig == 0){
 378           clusters[3][clno]=ncl[i];
 379         }else{
 380           clusters[3][clno]=cells[j];
 381         }
 382       }
 383     }
 384   }
 385
 386 }
 387
 388
 389 void AliPMDClustering::gaussfit(int ncell, int nclust, double &x, double &y ,double &z, double &xc, double &yc, double &zc, double &rc)
 390 {
 391   int i, j, i1, i2, jmax, novar, idd, jj;
 392   double xx[4500], yy[4500], zz[4500], xxc[4500], yyc[4500];
 393   double a[4500], b[4500], c[4500], d[4500], ha[4500], hb[4500];
 394   double hc[4500], hd[4500], zzc[4500], rrc[4500];
 395   int neib[4500][50];
 396   double sum, dx, dy, str, str1, aint, sum1, rr, dum;
 397   double x1, x2, y1, y2;
 398   str=0.;
 399   str1=0.;
 400   rr=0.3;
 401   novar=0;
 402
 403   j = 0;  // Just put not to see the compiler warning, BKN
 404
 405
 406   for(i=0; i<=ncell; i++){
 407     xx[i]=*(&x+i);
 408     yy[i]=*(&y+i);
 409     zz[i]=*(&z+i);
 410     str=str+zz[i];
 411   }
 412   for(i=0; i<=nclust; i++){
 413     xxc[i]=*(&xc+i);
 414     yyc[i]=*(&yc+i);
 415     zzc[i]=*(&zc+i);
 416     str1=str1+zzc[i];
 417     rrc[i]=0.5;
 418
 419   }
 420   for(i=0; i<=nclust; i++){
 421     zzc[i]=str/str1*zzc[i];
 422     ha[i]=xxc[i];
 423     hb[i]=yyc[i];
 424     hc[i]=zzc[i];
 425     hd[i]=rrc[i];
 426     x1=xxc[i];
 427     y1=yyc[i];
 428   }
 429   for(i=0; i<=ncell; i++){
 430     idd=0;
 431     x1=xx[i];
 432     y1=yy[i];
 433     for(j=0; j<=nclust; j++){
 434       x2=xxc[j];
 435       y2=yyc[j];
 436       if(Dist(x1,y1,x2,y2) <= 3.){ idd=idd+1; neib[i][idd]=j; }
 437     }
 438
 439     neib[i][0]=idd;
 440   }
 441   sum=0.;
 442   for(i1=0; i1<=ncell; i1++){
 443     aint=0.;
 444     idd=neib[i1][0];
 445     for(i2=1; i2<=idd; i2++){
 446       jj=neib[i1][i2];
 447       dx=xx[i1]-xxc[jj];
 448       dy=yy[i1]-yyc[jj];
 449       dum=rrc[j]*rrc[jj]+rr*rr;
 450       aint=aint+exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
 451     }
 452     sum=sum+(aint-zz[i1])*(aint-zz[i1])/str;
 453   }
 454   jmax=nclust*1000;
 455   if(nclust > 20)jmax=20000;
 456   for(j=0; j<jmax; j++){
 457     str1=0.;
 458     for(i=0; i<=nclust; i++){
 459       a[i]=xxc[i]+0.6*(ranmar()-0.5);
 460       b[i]=yyc[i]+0.6*(ranmar()-0.5);
 461       c[i]=zzc[i]*(1.+(ranmar()-0.5)*0.2);
 462       str1=str1+zzc[i];
 463       d[i]=rrc[i]*(1.+(ranmar()-0.5)*0.1);
 464       if(d[i] < 0.25)d[i]=0.25;
 465     }
 466     for(i=0; i<=nclust; i++){ c[i]=c[i]*str/str1; }
 467     sum1=0.;
 468     for(i1=0; i1<=ncell; i1++){
 469       aint=0.;
 470       idd=neib[i1][0];
 471       for(i2=1; i2<=idd; i2++){
 472         jj=neib[i1][i2];
 473         dx=xx[i1]-a[jj];
 474         dy=yy[i1]-b[jj];
 475         dum=d[jj]*d[jj]+rr*rr;
 476         aint=aint+exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
 477       }
 478       sum1=sum1+(aint-zz[i1])*(aint-zz[i1])/str;
 479     }
 480
 481     if(sum1 < sum){
 482       for(i2=0; i2<=nclust; i2++){
 483         xxc[i2]=a[i2];
 484         yyc[i2]=b[i2];
 485         zzc[i2]=c[i2];
 486         rrc[i2]=d[i2];
 487         sum=sum1;
 488
 489       }
 490     }
 491   }
 492   for(j=0; j<=nclust; j++){
 493     *(&xc+j)=xxc[j];
 494     *(&yc+j)=yyc[j];
 495     *(&zc+j)=zzc[j];
 496     *(&rc+j)=rrc[j];
 497   }
 498 }
 499
 500
 501 double AliPMDClustering::Dist(double x1, double y1, double x2, double y2)
 502 {
 503   return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 504 }
 505
 506
 507 int AliPMDClustering::crclust(double ave, double cutoff, int nmx1, int idet)
 508 {
 509   int i,j,k,id1,id2,icl, numcell, clust[2][5000];
 510   int jd1,jd2, icell, cellcount;
 511   static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 512   // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
 513   // cell. There are six neighbours.
 514   // cellcount --- total number of cells having nonzero ener dep
 515   // numcell --- number of cells in a given supercluster
 516   //ofstream ofl0("cells_loc",ios::out);
 517   // initialize infocl[2][ndimx][ndimy]
 518   for (j=0; j < 72; j++){
 519     for(k=0; k < 72; k++){
 520       infocl[0][j][k] = 0;
 521       infocl[1][j][k] = 0;
 522     }
 523   }
 524   for(i=0; i < nmx; i++){
 525     infcl[0][i] = -1;
 526     id1=iord[0][i];
 527     id2=iord[1][i];
 528     if(d[id1][id2] <= cutoff){infocl[0][id1][id2]=-1;}
 529   }
 530   // ---------------------------------------------------------------
 531   // crude clustering begins. Start with cell having largest adc
 532   // count and loop over the cells in descending order of adc count
 533   // ---------------------------------------------------------------
 534   icl=-1;
 535   cellcount=-1;
 536   for(icell=0; icell <= nmx1; icell++){
 537     id1=iord[0][icell];
 538     id2=iord[1][icell];
 539     if(infocl[0][id1][id2] == 0 ){
 540       // ---------------------------------------------------------------
 541       // icl -- cluster #, numcell -- # of cells in it, clust -- stores
 542       // coordinates of the cells in a cluster, infocl[0][i1][i2] is 1 for
 543       // primary and 2 for secondary cells,
 544       // infocl[1][i1][i2] stores cluster #
 545       // ---------------------------------------------------------------
 546       icl=icl+1;
 547       numcell=0;
 548       cellcount=cellcount+1;
 549       infocl[0][id1][id2]=1;
 550       infocl[1][id1][id2]=icl;
 551       infcl[0][cellcount]=icl;
 552       infcl[1][cellcount]=id1;
 553       infcl[2][cellcount]=id2;
 554       clust[0][numcell]=id1;
 555       clust[1][numcell]=id2;
 556       for(i=1; i<5000; i++)clust[0][i]=0;
 557       // ---------------------------------------------------------------
 558       // check for adc count in neib. cells. If ne 0 put it in this clust
 559       // ---------------------------------------------------------------
 560       for(i=0; i<6; i++){
 561         jd1=id1+neibx[i];
 562         jd2=id2+neiby[i];
 563         if( (jd1 >= 0 && jd1 < 72) && (jd2 >= 0 && jd2 < 72) &&
 564             infocl[0][jd1][jd2] == 0){
 565           numcell=numcell+1;
 566           infocl[0][jd1][jd2]=2;
 567           infocl[1][jd1][jd2]=icl;
 568           clust[0][numcell]=jd1;
 569           clust[1][numcell]=jd2;
 570           cellcount=cellcount+1;
 571           infcl[0][cellcount]=icl;
 572           infcl[1][cellcount]=jd1;
 573           infcl[2][cellcount]=jd2;
 574         }
 575       }
 576       // ---------------------------------------------------------------
 577       // check adc count for neighbour's neighbours recursively and
 578       // if nonzero, add these to the cluster.
 579       // ---------------------------------------------------------------
 580       for(i=1;i < 5000;i++){
 581         if(clust[0][i] != 0){
 582           id1=clust[0][i];
 583           id2=clust[1][i];
 584           for(j=0; j<6 ; j++){
 585             jd1=id1+neibx[j];
 586             jd2=id2+neiby[j];
 587             if( (jd1 >= 0 && jd1 < 72) && (jd2 >= 0 && jd2 < 72) &&
 588                 infocl[0][jd1][jd2] == 0 ){
 589               infocl[0][jd1][jd2]=2;
 590               infocl[1][jd1][jd2]=icl;
 591               numcell=numcell+1;
 592               clust[0][numcell]=jd1;
 593               clust[1][numcell]=jd2;
 594               cellcount=cellcount+1;
 595               infcl[0][cellcount]=icl;
 596               infcl[1][cellcount]=jd1;
 597               infcl[2][cellcount]=jd2;
 598             }
 599           }
 600         }
 601       }
 602     }
 603   }
 604   //  for(icell=0; icell<=cellcount; icell++){
 605   //    ofl0 << infcl[0][icell] << " " << infcl[1][icell] << " " <<
 606   //      infcl[2][icell] << endl;
 607   //  }
 608   return cellcount;
 609 }
 610
 611 double AliPMDClustering::ranmar()
 612 {
 613   /*                                   C==========================C*/
 614   /*===================================C==========================*/
 615   /*  Universal random number generator proposed by Marsaglia and Zaman
 616       in report FSU-SCRI-87-50 */
 617
 618   //  clock_t start;
 619   int ii, jj;
 620   static int i=96, j=32, itest=0, i1, i2, i3, i4, i5;
 621   static double u[97], c, cd, cm, s, t;
 622   static double uni;
 623   int count1,count2,idum;
 624   /*    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  */
 625   if (itest == 0) {
 626     //*******************************************************
 627     // following three lines if the seed to be provided by computer
 628     // start = time(NULL);
 629     // ii=start;
 630     // jj=start;
 631     //*******************************************************
 632     //following two lines for fixed seed ( during testing only. Else
 633     //use preceeing three lines
 634     ii=8263;
 635     jj=5726;
 636     if(ii > 31328 ) ii = ii - ( ii / 31328 ) * 31328;
 637     if(jj > 30081 ) jj = jj - ( jj / 30081 ) * 30081;
 638     itest=itest+1;
 639     if((( ii > 0 ) &&  ( ii <= 31328 )) && (( jj > 0 ) &&
 640                                             ( jj <= 30081 ))){
 641       i1=ii/177+2; i2=ii-(i1-2)*177+2; i3=jj/169+1; i4=jj-(i3-1)*169;
 642       i4 = jj - (i3-1)*169;
 643       count1=0;
 644       while ( count1 < 97 ){
 645         s=0.;
 646         t=0.5;
 647         count2=0;
 648         while( count2 < 24 ){
 649           idum=i1*i2/179;
 650           idum=( i1*i2 - (i1*i2/179)*179 ) * i3;
 651           i5=idum-(idum/179)*179;
 652           i1=i2; i2=i3; i3=i5; idum=53*i4+1; i4=idum-(idum/169)*169;
 653           if( i4*i5-((i4*i5)/64)*64 >= 32 ) s=s+t;
 654           t=0.5*t;
 655           count2=count2+1;
 656         }
 657         u[count1] = s;
 658         count1 = count1 +1;
 659       }
 660       c = 362436./16777216.;  cd = 7654321./16777216.;
 661       cm = 16777213./16777216.;
 662     }
 663     else{
 664       cout << " wrong initialization " << endl;
 665     }
 666   }
 667   else{
 668     uni = u[i] - u[j]; if( uni < 0.) uni = uni + 1; u[i] = uni;
 669     i = i -1;
 670     if( i < 0 ) i = 96; j = j - 1; if ( j < 0 ) j = 96; c = c - cd;
 671     if( c < 0. ) c = c+cm; uni = uni-c ; if( uni < 0. )uni = uni+1.;
 672     //    return uni;
 673   }
 674   return uni;
 675
 676 }
 677
 678 void AliPMDClustering::ConvertL2G(int smnumber, double xcell, double ycell, double &xpos, double &ypos)
 679 {
 680   float xreal = -999., yreal = -999.;
 681   float cell_rad=0.25, celldia_x=0.5, celldia_y=0.4330127;
 682   float xcon, ycon;
 683   float xoff1, xoff2, yoff=0.2886751, yoff3;
 684   float xhex1 = -27.09375, yhex1 = -15.652584;
 685   float xhex2 =  27.09375, yhex2 = -15.652584;
 686   float xhex3 = 0.0, yhex3 = 31.285168;
 687
 688
 689   double xcorner[27] =
 690   {
 691     9.435395, 45.560394,  81.685394, -8.627106,
 692     27.497894, 63.622894, -26.689606,  9.435394,
 693     45.560394, 9.435344,  -8.627106, -26.689556,
 694     45.560345, 27.497894, 9.435445,   81.685341,
 695     63.622894, 45.560444, -18.870789, -36.933388,
 696     -54.995991, -36.933189, -54.995789, -73.058388,
 697     -54.995586, -73.058189, -91.120789
 698   };
 699
 700   double ycorner[27] =
 701   {
 702     -16.342583, -16.34258, -16.34258, -47.627750, -47.627750,
 703     -47.627750, -78.912918, -78.912918, -78.912918, 16.342611,
 704     47.627808, 78.913002, 16.342554, 47.627750, 78.912949,
 705     16.342495, 47.627693, 78.912888, -0.000116, -31.285227,
 706     -62.570335, 31.285110, 0.000000, -31.285110, 62.570335,
 707     31.285227, 0.000116
 708   };
 709
 710   if (smnumber<=8)
 711     {
 712       xcon  = xcorner[smnumber]+xhex1;
 713       ycon  = ycorner[smnumber]+yhex1;
 714       xoff1 = celldia_x+(ycell-1)*cell_rad;
 715       xreal = xcon+xoff1+celldia_x*(xcell-1);
 716       yreal = ycon+yoff+celldia_y*(ycell-1);
 717     }
 718
 719   if (smnumber>8 && smnumber<=17)
 720     {
 721       xcon  = xcorner[smnumber]+xhex2;
 722       ycon  = ycorner[smnumber]+yhex2;
 723       xoff2 = celldia_x+(xcell-1)*cell_rad;
 724       xreal = xcon-(xoff2+celldia_x*(ycell-1));
 725       yreal = ycon+yoff+celldia_y*(xcell-1);
 726     }
 727
 728   if (smnumber>17)
 729     {
 730       xcon  = xcorner[smnumber]+xhex3;
 731       ycon  = ycorner[smnumber]+yhex3;
 732       yoff3 = celldia_x * 0.8660254  + cell_rad * 0.57735027;
 733       xreal = xcon+(ycell-xcell)*cell_rad;
 734       yreal = ycon-(yoff3+(xcell+ycell-2)*celldia_y);
 735     }
 736
 737   xpos = xreal;
 738   ypos = yreal;
 739 }
 740
 741 void AliPMDClustering::cell_pos(Int_t isup, Int_t j, int k, Float_t &xp, Float_t &yp){
 742
 743   /*
 744     This converts PMD cluster or CELL coordinates
 745     to Global coordinates.
 746     Written by Prof. S.C. Phatak
 747   */
 748
 749   Int_t i;
 750   Float_t celldia = 0.5;
 751   const Float_t pi = 3.14159;
 752   const double sqrth=0.8660254;  // sqrth = sqrt(3.)/2.
 753   /*
 754     isup --> supermodule no ( 0 - 26 )
 755     idet --> detector ( pmd or cpv : not required now )
 756     j --> xpad ( goes from 1 to 72 )
 757     k --> ypad ( goes from 1 to 72 )
 758     xp --> global x coordinate
 759     yp --> global y coordinate
 760
 761     (xp0,yp0) corner positions of all supermodules in global
 762     coordinate system. That is the origin
 763     of the local ( supermodule ) coordinate system.
 764 */
 765
 766   Float_t xp0[27] =
 767   {
 768     -17.9084, 18.2166, 54.3416, -35.9709, 0.154144,
 769     36.2791, -54.0334, -17.9084, 18.2166, 36.7791,
 770     18.7166, 0.654194, 72.9041, 54.8416, 36.7792,
 771     109.029, 90.9666, 72.9042, -18.8708, -36.9334,
 772     -54.996, -36.9332, -54.9958, -73.0584, -54.9956,
 773     -73.0582, -91.1208
 774   };
 775
 776   Float_t yp0[27] =
 777   {
 778     -32.1395, -32.1395, -32.1395, -63.4247, -63.4247,
 779     -63.4247, -94.7098, -94.7098, -94.7098, 0.545689,
 780     31.8309, 63.1161, 0.545632, 31.8308, 63.116,
 781     0.545573, 31.8308, 63.116, 31.5737, 0.288616,
 782     -30.9965, 62.859, 31.5738, 0.288733, 94.1442,
 783     62.8591, 31.574
 784   };
 785
 786   /*
 787      angles of rotation for three sets of supermodules
 788      The angle is same for first nine, next nine and last nine
 789      supermodules
 790   */
 791
 792   Float_t th[3] = {0., -2.*pi/3., 2.*pi/3.};
 793   Float_t xr, yr, xinit, yinit, cs, sn;
 794
 795   /*
 796      xinit and yinit are coordinates of the cell in local coordinate system
 797   */
 798
 799   xinit = (j)*celldia+(k)/2.*celldia;
 800   yinit = sqrth*(k)/2.;
 801   i=isup/9;
 802   cs=cos(th[i]);
 803   sn=sin(th[i]);
 804   //
 805   // rotate first
 806   //
 807   xr=cs*xinit+sn*yinit;
 808   yr=-sn*xinit+cs*yinit;
 809   //
 810   // then translate
 811   //
 812   xp=xr+xp0[isup];
 813   yp=yr+yp0[isup];
 814
 815 }
 816 void AliPMDClustering::SetMessage(Int_t imessage)
 817 {
 818   fMessage = imessage;
 819 }