HLT/CALO/AliHLTCaloClusterizer.cxx

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