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