HLT/PHOS/AliHLTPHOSClusterizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Authors: Øystein Djuvsland <oysteind@ift.uib.no>                       *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15
  16 /** @file   AliHLTPHOSClusterizer.cxx
  17     @author Øystein Djuvsland
  18     @date
  19     @brief  A temporary clusterizer for PHOS
  20 */
  21
  22
  23
  24 //#include "AliHLTPHOSDefinitions.h"
  25 #include "AliHLTPHOSClusterizer.h"
  26 //#include "AliHLTPHOSCommonDefs.h"
  27 #include "TVector3.h"
  28 #include "TMath.h"
  29 #include "AliHLTPHOSRcuCellEnergyDataStruct.h"
  30 #include "AliHLTPHOSRecPointListDataStruct.h"
  31 #include "AliHLTPHOSValidCellDataStruct.h"
  32 #include "AliHLTPHOSRecPointDataStruct.h"
  33 #include "AliHLTPHOSClusterDataStruct.h"
  34 //#include "AliHLTPHOSConstants.h"
  35
  36 //using namespace PhosHLTConst;
  37
  38
  39
  40 ClassImp(AliHLTPHOSClusterizer);
  41
  42 /**
  43 * Main constructor
  44 **/
  45 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():AliHLTPHOSBase(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0),
  46                                                fHighGainFactor(0.005), fLowGainFactor(0.08),
  47                                                fArraySize(3), fMultiplicity(fArraySize*fArraySize)
  48 //AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():AliHLTPHOSProcessor(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0),
  49 //                                             fHighGainFactor(0.005), fLowGainFactor(0.08),
  50 //                                             fArraySize(3), fMultiplicity(fArraySize*fArraySize)
  51 {
  52   //See header file for documentation
  53
  54 }//end
  55
  56 // PTH AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):AliHLTPHOSBase(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0),
  57 //                                                                          fHighGainFactor(0.005), fLowGainFactor(0.08),
  58 //                                                                          fArraySize(3), fMultiplicity(fArraySize*fArraySize)
  59
  60 //{
  61   //Copy constructor, not implemented
  62 //}//end
  63
  64 AliHLTPHOSClusterizer:: ~AliHLTPHOSClusterizer()
  65 {
  66   //Destructor
  67 }
  68
  69 /**
  70 * Building a 2D array of cell energies of the PHOS detector
  71 * @param cellData object containing the cell energies from one event
  72 * @param recPointList list to be filled with coordinates of local maxima in the detector
  73 **/
  74
  75 int
  76 AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData,
  77                                             AliHLTPHOSRecPointListDataStruct* recPointList)
  78 {
  79   //BuildCellEnergyArray
  80
  81   Int_t x = 0;
  82   Int_t z = 0;
  83   Int_t gain = 0;
  84   Int_t xMod = 0;
  85   Int_t zMod = 0;
  86   Int_t index = 0;
  87   Double_t energyCount = 0;
  88
  89   for(Int_t cell = 0; cell < cellData->fCnt; cell++)
  90     {
  91
  92       z = (cellData->fValidData[cell]).fZ;
  93       x = (cellData->fValidData[cell]).fX;
  94       gain = (cellData->fValidData[cell]).fGain;
  95
  96       zMod = z + (cellData->fRcuZ)*N_ZROWS_RCU;
  97       xMod = x + (cellData->fRcuX)*N_XCOLUMNS_RCU;
  98
  99       energyCount = (cellData->fValidData[cell]).fEnergy;
 100
 101       if(gain == 0 && energyCount < 1023)
 102         {
 103           fEnergyArray[xMod][zMod] = fHighGainFactor * energyCount;
 104
 105           if(fEnergyArray[xMod][zMod] > fClusterThreshold)
 106             {
 107               recPointList[index].fZ = zMod;
 108               recPointList[index].fX = xMod;
 109
 110               for(Int_t j = 0; j < index; j++)
 111                 {
 112                   if(recPointList[j].fZ == zMod && recPointList[j].fX == xMod)
 113                     recPointList[j].fZ = -1;
 114                 }
 115               index++;
 116             }
 117
 118         }
 119       else if(gain == 1 && fEnergyArray[xMod][zMod] == 0)
 120         {
 121           fEnergyArray[xMod][zMod] = fLowGainFactor * energyCount;
 122            if(fEnergyArray[xMod][zMod] > fClusterThreshold)
 123             {
 124               recPointList[index].fZ = zMod;
 125               recPointList[index].fX = xMod;
 126               recPointList[index].fModule = cellData->fModuleID;
 127               index++;
 128             }
 129         }
 130
 131     }
 132
 133   fPHOSModule = cellData->fModuleID;
 134
 135   return index;
 136 }//end BuildCellEnergyArray
 137
 138
 139
 140 /**
 141 * Creating an array of rec points
 142 * @param recPointStructArrayPtr array to store the rec points
 143 * @param list list of rec points
 144 * @param nPoints number of points
 145 **/
 146 int
 147 AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr,
 148                                                  AliHLTPHOSRecPointListDataStruct* list,
 149                                                  int nPoints)
 150
 151 {
 152   //CreateRecPointStructArray
 153
 154   Int_t flag = 0;
 155   Int_t edgeFlagRows = 0;
 156   Int_t edgeFlagCols = 0;
 157   Int_t k = 0;
 158   Int_t nRecPoints = 0;
 159   Int_t z = 0;
 160   Int_t x = 0;
 161   Int_t i = 0;
 162   Int_t j = 0;
 163   Int_t a = 0;
 164
 165   Float_t* energiesList = NULL;
 166
 167   for(Int_t point = 0; point < nPoints; point++)
 168     {
 169       z = list[point].fZ;
 170       x = list[point].fX;
 171       if(z == -1) continue;
 172
 173       if((z-fArraySize/2) < 0/*= - N_ZROWS_MOD/2*/ || (z+fArraySize/2) >= N_ZROWS_MOD/*/2*/)
 174         {
 175           edgeFlagRows = 1;
 176           continue;
 177         }
 178       if((x-fArraySize/2) < 0/*= - N_COLUMNS_MOD/2*/ || (x+fArraySize/2) >= N_XCOLUMNS_MOD)
 179         {
 180           edgeFlagCols = 1;
 181           continue;
 182         }
 183
 184
 185       if(!flag) energiesList = new Float_t[fMultiplicity];
 186
 187       flag = 0;
 188       k = 0;
 189
 190       for(i = -fArraySize/2; i <= fArraySize/2; i++)
 191         {
 192           if(flag) break;
 193           for(j = -fArraySize/2; j <= fArraySize/2; j++)
 194             {
 195
 196               if(fEnergyArray[x+i][z+j] > fEnergyArray[x][z] && abs(i) < 2 && abs(j) < 2)
 197                 {
 198                   flag = 1;
 199                   break;
 200                 }
 201
 202               energiesList[k] = fEnergyArray[x+i][z+j];
 203               k++;
 204             }
 205         }
 206
 207
 208       if(!flag && k)
 209         {
 210           for(a = 0; a < k; a++)
 211             {
 212               (recPointStructArrayPtr[nRecPoints].fEnergiesListPtr)[a] = energiesList[a];
 213             }
 214           recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[0] = x;
 215           recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[1] = z;
 216           recPointStructArrayPtr[nRecPoints].fX = x;
 217           recPointStructArrayPtr[nRecPoints].fZ = z;
 218           // recPointStructArrayPtr[nRecPoints].fMultiplicity = fMultiplicity;
 219           recPointStructArrayPtr[nRecPoints].fPHOSModule = list[point].fModule;
 220           nRecPoints++;
 221         }
 222     }
 223   if(energiesList)
 224     {
 225       delete [] energiesList;
 226       energiesList = NULL;
 227     }
 228
 229   return nRecPoints;
 230
 231 }//end CreateRecPointStructArray
 232
 233
 234 /**
 235 * Calculating the center of gravity of a rec point
 236 * Not working well at this point!
 237 * @param recPointPtr pointer to the rec point
 238 **/
 239 int
 240 AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
 241 {
 242   //CalculateCenterOfGravity
 243   //Copied from offline code, and modified
 244
 245   Float_t xt = 0;
 246   Float_t zt = 0;
 247   Float_t xi = 0;
 248   Float_t zj = 0;
 249   Float_t w = 0;
 250   Float_t w0 = 4.5;
 251   Float_t wtot = 0;
 252
 253   Int_t x = recPointPtr->fCoordinatesPtr[0];
 254   Int_t z = recPointPtr->fCoordinatesPtr[1];
 255
 256   for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
 257     {
 258       for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
 259         {
 260           xi = x + i;
 261           zj = z + j;
 262           w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
 263           //printf("log in cog: %f\n", log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]));
 264           xt += xi * w ;
 265           zt += zj * w ;
 266           wtot += w ;
 267         }
 268     }
 269   xt /= wtot ;
 270   zt /= wtot ;
 271
 272   //printf("wtot cog: %f\n", wtot);
 273
 274   recPointPtr->fX = xt;
 275   recPointPtr->fZ = zt;
 276
 277   return 0;
 278
 279 }//end CalculateCenterOfGravity
 280
 281 int
 282 AliHLTPHOSClusterizer::CalculateMoments(AliHLTPHOSRecPointDataStruct* recPointPtr, Bool_t axisOnly)
 283 {
 284   //Copied from offline code, and modified
 285
 286
 287   // Calculate the shower moments in the eigen reference system
 288   // M2x, M2z, M3x, M4z
 289   // Calculate the angle between the shower position vector and the eigen vector
 290
 291   Double_t wtot = 0. ;
 292   Double_t x    = 0.;
 293   Double_t z    = 0.;
 294   Double_t dxx  = 0.;
 295   Double_t dzz  = 0.;
 296   Double_t dxz  = 0.;
 297   Double_t lambda0=0, lambda1=0;
 298   Float_t logWeight = 4.5;
 299
 300   //Int_t iDigit;
 301
 302   // 1) Find covariance matrix elements:
 303   //    || dxx dxz ||
 304   //    || dxz dzz ||
 305
 306   Int_t xi;
 307   Int_t zj;
 308   Double_t w;
 309
 310   Int_t xc = recPointPtr->fCoordinatesPtr[0];
 311   Int_t zc = recPointPtr->fCoordinatesPtr[1];
 312
 313   for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
 314     {
 315       xi = xc + i;
 316       for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
 317         {
 318           zj = zc + j;
 319           if (fEnergyArray[xi][zj] > 0)
 320             {
 321
 322               //              w     = TMath::Max(0.,logWeight + log( fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
 323
 324               //printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
 325               x    += w * xi ;
 326               z    += w * zj ;
 327               dxx  += w * xi * xi ;
 328               dzz  += w * zj * zj ;
 329               dxz  += w * xi * zj ;
 330               wtot += w ;
 331             }
 332         }
 333     }
 334   //printf("wtot: %f\n", wtot);
 335   if (wtot>0)
 336     {
 337       x   /= wtot ;
 338       z   /= wtot ;
 339       dxx /= wtot ;
 340       dzz /= wtot ;
 341       dxz /= wtot ;
 342       dxx -= x * x ;
 343       dzz -= z * z ;
 344       dxz -= x * z ;
 345
 346   // 2) Find covariance matrix eigen values lambda0 and lambda1
 347
 348       lambda0 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 349       lambda1 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 350     }
 351
 352   recPointPtr->fM2x   = lambda0;
 353   recPointPtr->fM2z   = lambda1;
 354   //  printf("lambda0 = %f -- lambda1 = %f\n", lambda0, lambda1);
 355   // 3) Find covariance matrix eigen vectors e0 and e1
 356   if(!axisOnly)
 357     {
 358       TVector2 e0,e1;
 359       if (dxz != 0)
 360         e0.Set(1.,(lambda0-dxx)/dxz);
 361       else
 362         e0.Set(0.,1.);
 363
 364       e0 = e0.Unit();
 365       e1.Set(-e0.Y(),e0.X());
 366
 367       // 4) Rotate cluster tensor from (x,z) to (e0,e1) system
 368       //    and calculate moments M3x and M4z
 369
 370       Float_t cosPhi = e0.X();
 371       Float_t sinPhi = e0.Y();
 372
 373       Float_t xiPHOS ;
 374       Float_t zjPHOS ;
 375       Double_t dx3, dz3, dz4;
 376       wtot = 0.;
 377       x    = 0.;
 378       z    = 0.;
 379       dxx  = 0.;
 380       dzz  = 0.;
 381       dxz  = 0.;
 382       dx3  = 0.;
 383       dz3  = 0.;
 384       dz4  = 0.;
 385       for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
 386         {
 387           for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
 388             {
 389               xi = xc + i;
 390               zj = zc + j;
 391               xiPHOS = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
 392               zjPHOS = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
 393               //              xi = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
 394               //              zj = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
 395               if (fEnergyArray[xi][zj] > 0)
 396                 {
 397                   w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
 398                   //  printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
 399
 400                   x    += w * xi ;
 401                   z    += w * zj ;
 402                   dxx  += w * xi * xi ;
 403                   dzz  += w * zj * zj ;
 404                   dxz  += w * xi * zj ;
 405                   dx3  += w * xi * xi * xi;
 406                   dz3  += w * zj * zj * zj ;
 407                   dz4  += w * zj * zj * zj * zj ;
 408                   wtot += w ;
 409                 }
 410             }
 411         }
 412       if (wtot>0)
 413         {
 414           x   /= wtot ;
 415           z   /= wtot ;
 416           dxx /= wtot ;
 417           dzz /= wtot ;
 418           dxz /= wtot ;
 419           dx3 /= wtot ;
 420           dz3 /= wtot ;
 421           dz4 /= wtot ;
 422           dx3 += -3*dxx*x + 2*x*x*x;
 423           dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z;
 424           dxx -= x * x ;
 425           dzz -= z * z ;
 426           dxz -= x * z ;
 427         }
 428
 429       // 5) Find an angle between cluster center vector and eigen vector e0
 430
 431       Float_t phi = 0;
 432       if ( (x*x+z*z) > 0 )
 433         phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z));
 434
 435       recPointPtr->fM3x   = dx3;
 436       recPointPtr->fM4z   = dz4;
 437       recPointPtr->fPhixe = phi;
 438       //      printf("%f\n",phi);
 439     }
 440   return 0;
 441 }
 442
 443
 444 /**
 445 * Create a simpler data struct for a rec point, containing
 446 * only coordinates, energy and timing
 447 * @param recPointPtr pointer to the rec point
 448 * @clusterStructPtr pointer to the cluster struct
 449 **/
 450 int
 451 AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
 452 {
 453   //ClusterizeStruct
 454
 455   Float_t clusterEnergy = 0;
 456   Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
 457
 458   for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
 459     {
 460       clusterEnergy += energiesListPtr[i];
 461     }
 462
 463   clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
 464   clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
 465   clusterStructPtr->fClusterEnergy = clusterEnergy;
 466   clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;
 467
 468   return 0;
 469
 470 }//end ClusterizeStruct
 471
 472
 473
 474
 475 /**
 476 * Resets the cell energy array
 477 **/
 478 int
 479 AliHLTPHOSClusterizer::ResetCellEnergyArray()
 480
 481 {
 482   //ResetCellEnergyArray
 483
 484   for(Int_t x = 0; x < N_ZROWS_MOD; x++)
 485     {
 486       for(Int_t z = 0; z < N_XCOLUMNS_MOD; z++)
 487         {
 488           fEnergyArray[x][z] = 0;
 489         }
 490     }
 491
 492   return 0;
 493
 494 }//end ResetCellEnergyArray
 495