HLT/CALO/AliHLTCaloClusterizer.cxx

   1 // $Id$
   2
   3 /**************************************************************************
   4  * This file is property of and copyright by the ALICE HLT Project        *
   5  * All rights reserved.                                                   *
   6  *                                                                        *
   7  * Primary Authors: Oystein Djuvsland                                     *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 /**
  19  * @file   AliHLTCaloClusterizer.cxx
  20  * @author Oystein Djuvsland
  21  * @date
  22  * @brief  Clusterizer for PHOS HLT
  23  */
  24
  25 // see header file for class documentation
  26 // or
  27 // refer to README to build package
  28 // or
  29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  30
  31 #include "AliHLTCaloClusterizer.h"
  32 #include "AliHLTLogging.h"
  33 #include "TMath.h"
  34 #include "AliHLTCaloRecPointDataStruct.h"
  35 #include "AliHLTCaloDigitDataStruct.h"
  36 #include "AliHLTCaloDigitContainerDataStruct.h"
  37 #include "AliHLTCaloConstantsHandler.h"
  38
  39 ClassImp(AliHLTCaloClusterizer);
  40
  41 AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):
  42         AliHLTCaloConstantsHandler(det),
  43         fCompareFunction(CompareDigitsByPosition),
  44         fRecPointArray(0),
  45         fRecPointDataPtr(0),
  46         fFirstRecPointPtr(0),
  47         fArraySize(0),
  48         fAvailableSize(0),
  49         fUsedSize(0),
  50         fNRecPoints(0),
  51         fDigitIndexPtr(0),
  52         fEmcClusteringThreshold(0),
  53         fEmcMinEnergyThreshold(0),
  54         fEmcTimeGate(0),
  55         fDigitsInCluster(0),
  56         fDigitsPointerArray(0),
  57         fDigitContainerPtr(0),
  58         fMaxDigitIndexDiff(0),
  59         fNDigits(0),
  60         fSortedByPosition(false),
  61         fSortedByEnergy(false),
  62         fSortDigits(false),
  63         fIsEMCAL(false),
  64         fBuffer(0)
  65
  66 {
  67     //See header file for documentation
  68     //fEmcClusteringThreshold = 0.2;
  69     //fEmcMinEnergyThreshold = 0.03;
  70
  71     fEmcClusteringThreshold = 0.1;
  72     fEmcMinEnergyThreshold = 0.01;
  73     fEmcTimeGate = 1.e-6 ;
  74
  75     fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();
  76
  77
  78     fArraySize = 10;
  79     fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];
  80
  81     fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;
  82     fBuffer = new UChar_t[fAvailableSize]; //FR
  83
  84     fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);
  85     fRecPointDataPtr = fFirstRecPointPtr;
  86
  87 }//end
  88
  89 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()
  90 {
  91     //See header file for documentation
  92   delete [] fBuffer; //FR
  93   fBuffer = NULL;    //FR
  94   fAvailableSize = 0;  //FR
  95
  96   delete [] fRecPointArray;
  97 }
  98
  99 void
 100 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)
 101 {
 102     // See header file for documentation
 103     fRecPointDataPtr = recPointDataPtr;
 104 }
 105
 106 Int_t
 107 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)
 108 {
 109     //see header file for documentation
 110     Int_t nRecPoints = 0;
 111     fNRecPoints = 0;
 112     fUsedSize = 0;
 113     fNDigits = nDigits;
 114     fRecPointDataPtr = fFirstRecPointPtr;
 115
 116     // Sort our digits
 117     SortDigits();
 118
 119     //Clusterization starts
 120     for (Int_t i = 0; i < nDigits; i++)
 121     {
 122         fDigitsInCluster = 0;
 123
 124          HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);
 125
 126         if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)
 127         {
 128            // Since we have sorted by energy the next digit will have even lower energy, so we return
 129            return fNRecPoints;
 130         }
 131
 132         if(fDigitsPointerArray[i]->fAssociatedCluster != -1)
 133         {
 134            // The digit is added to a previous cluster, continue
 135            continue;
 136         }
 137
 138         CheckArray();
 139         CheckBuffer();
 140
 141         // First digit is placed at the fDigits member variable in the recpoint
 142         fDigitIndexPtr = &(fRecPointDataPtr->fDigits);
 143
 144         fRecPointDataPtr->fAmp = 0;
 145         fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;
 146
 147         // Assigning the digit to this rec point
 148         fRecPointDataPtr->fDigits = i;
 149         fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);
 150
 151         // Incrementing the pointer to be ready for new entry
 152         fDigitIndexPtr++;
 153
 154         fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;
 155
 156
 157         //fDigitsPointerArray[i]->fEnergy = 0;
 158         fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;
 159
 160
 161         fDigitsInCluster++;
 162         nRecPoints++;
 163
 164         // Scanning for the neighbours
 165         if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)
 166         {
 167             return -1;
 168         }
 169
 170         //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);
 171
 172         fRecPointDataPtr->fMultiplicity = fDigitsInCluster;
 173         fRecPointArray[fNRecPoints] = fRecPointDataPtr;
 174
 175         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);
 176
 177         fNRecPoints++;
 178
 179     }//end of clusterization
 180
 181     return nRecPoints;
 182 }
 183
 184 Int_t
 185 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)
 186 {
 187     //see header file for documentation
 188
 189     // The following cuts can be used if we sort by posisiton. Not tested, but it should be fine...
 190     Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);
 191     Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));
 192
 193     // All digits for now
 194     max = fNDigits;
 195     min = 0;
 196
 197     for (Int_t j = min; j < max; j++)
 198     {
 199         if (fDigitsPointerArray[j]->fAssociatedCluster == -1 &&  fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)
 200         {
 201             if (j != index)
 202             {
 203                 if (AreNeighbours(fDigitsPointerArray[index],
 204                                   fDigitsPointerArray[j]))
 205                 {
 206                     // Check that the buffer is large enough for adding a digit (can be heavily improved wrt performance)
 207                     CheckBuffer();
 208
 209                     // Assigning index to digit
 210                     *fDigitIndexPtr = j;
 211                     fUsedSize += sizeof(Int_t);
 212
 213                     // Incrementing digit pointer to be ready for new entry
 214                     fDigitIndexPtr++;
 215
 216                     // Adding the digit energy to the rec point
 217                     fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;
 218
 219                     // Setting energy to 0
 220                     //fDigitsPointerArray[j]->fEnergy = 0;
 221
 222                     // Setting the associated cluster
 223                     fDigitsPointerArray[j]->fAssociatedCluster = fNRecPoints;
 224
 225                     HLTDebug("Added digit with index: %d, energy: %f, to associated cluster: %d", fDigitsPointerArray[j]->fID, fDigitsPointerArray[j]->fEnergy, fDigitsPointerArray[j]->fAssociatedCluster);
 226
 227                     fDigitsInCluster++;
 228
 229                     // Scan for neighbours of this digit
 230                     ScanForNeighbourDigits(j, recPoint);
 231                 }
 232             }
 233         }
 234     }
 235     return 0;
 236 }
 237
 238 Int_t
 239 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,
 240                                      AliHLTCaloDigitDataStruct* digit2)
 241 {
 242     //see header file for documentation
 243     if ( (digit1->fModule == digit2->fModule) || AreEdgeCells(digit1, digit2))
 244     {
 245         Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
 246         Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
 247
 248         // Common edge defines neighbour
 249         //if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))
 250           // Common edge and corner defines neighbour
 251         if (( coldiff <= 1   &&  rowdiff <= 1 ))
 252         {
 253             // Check also for time
 254             if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
 255             {
 256                 return 1;
 257             }
 258         }
 259     }
 260     return 0;
 261 }
 262
 263
 264
 265 Int_t AliHLTCaloClusterizer::CheckArray()
 266 {
 267     // See header file for class documentation
 268     if (fArraySize == fNRecPoints)
 269     {
 270         fArraySize *= 2;
 271         AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];
 272         memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));
 273         delete [] fRecPointArray;
 274         fRecPointArray = tmp;
 275     }
 276     return 0;
 277 }
 278
 279 Int_t AliHLTCaloClusterizer::CheckBuffer()
 280 {
 281     // See header file for class documentation
 282     if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))
 283     {
 284         Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);
 285         Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);
 286         UChar_t *tmp = new UChar_t[fAvailableSize*2];
 287
 288         if (tmp == NULL)
 289           {
 290             HLTError("Pointer error");
 291             return(-1);
 292           }
 293
 294         memcpy(tmp, fFirstRecPointPtr, fUsedSize);
 295         fAvailableSize *= 2;
 296         for (Int_t n = 0; n < fNRecPoints; n++)
 297         {
 298             fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));
 299         }
 300         delete [] fBuffer; //FR
 301         fBuffer = tmp; //FR
 302
 303         fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);
 304         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);
 305         fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);
 306         //fUsedSize = 0;
 307     }
 308     return 0;
 309 }
 310
 311 void AliHLTCaloClusterizer::SetSortDigitsByPosition()
 312 {
 313     // Sort the digit pointers by position
 314     fCompareFunction = &CompareDigitsByPosition;
 315     fSortDigits = true;
 316     fSortedByPosition = true;
 317 }
 318
 319 void AliHLTCaloClusterizer::SetSortDigitsByEnergy()
 320 {
 321     // See header file for class documentation
 322     fCompareFunction = &CompareDigitsByEnergy;
 323     fSortDigits = true;
 324     fSortedByEnergy = true;
 325 }
 326
 327 void AliHLTCaloClusterizer::SortDigits()
 328 {
 329     // See header file for class documentation
 330     if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);
 331 }
 332
 333 Int_t
 334 AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)
 335 {
 336     // See header file for documentation
 337     return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;
 338 }
 339
 340 Int_t
 341 AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)
 342 {
 343     // See header file for documentation
 344   if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;
 345   return 1;
 346 }
 347
 348 void AliHLTCaloClusterizer::SetDetector(TString det)
 349 {
 350   if(det.CompareTo("EMCAL"))
 351   {
 352     fIsEMCAL = true;
 353   }
 354   else
 355   {
 356     fIsEMCAL = false;
 357   }
 358 }
 359
 360 Bool_t AliHLTCaloClusterizer::AreEdgeCells(AliHLTCaloDigitDataStruct *digit0, AliHLTCaloDigitDataStruct *digit1)
 361 {
 362   if(fIsEMCAL)
 363   {
 364     Int_t modDiff = digit0->fModule - digit1->fModule;
 365     if(TMath::Abs(modDiff) > 1) return kFALSE;
 366     if(digit0->fModule > digit1->fModule && digit1->fModule%2 == 0)
 367     {
 368       if(digit0->fZ == 0 && digit1->fZ == (fCaloConstants->GetNZROWSMOD()-1))
 369       return kTRUE;
 370     }
 371     if(digit1->fModule > digit0->fModule && digit0->fModule%2 == 0)
 372     {
 373       if(digit1->fZ == 0 && digit0->fZ == (fCaloConstants->GetNZROWSMOD()-1))
 374       return kTRUE;
 375     }
 376   }
 377
 378   return false;
 379
 380 }