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         fBuffer(0)\r
  64 {\r
  65     //See header file for documentation\r
  66     //fEmcClusteringThreshold = 0.2;\r
  67     //fEmcMinEnergyThreshold = 0.03;\r
  68 \r
  69     fEmcClusteringThreshold = 0.1;\r
  70     fEmcMinEnergyThreshold = 0.01;\r
  71     fEmcTimeGate = 1.e-6 ;\r
  72 \r
  73     fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
  74 \r
  75 \r
  76     fArraySize = 10;\r
  77     fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
  78 \r
  79     fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
  80     fBuffer = new UChar_t[fAvailableSize]; //FR\r
  81     \r
  82     fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);\r
  83     fRecPointDataPtr = fFirstRecPointPtr;\r
  84 \r
  85 }//end\r
  86 \r
  87 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()\r
  88 {\r
  89     //See header file for documentation\r
  90   delete [] fBuffer; //FR\r
  91   fBuffer = NULL;    //FR\r
  92   fAvailableSize = 0;  //FR\r
  93 \r
  94   delete [] fRecPointArray;\r
  95 }\r
  96 \r
  97 void\r
  98 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
  99 {\r
 100     // See header file for documentation\r
 101     fRecPointDataPtr = recPointDataPtr;\r
 102 }\r
 103 \r
 104 Int_t\r
 105 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
 106 {\r
 107     //see header file for documentation\r
 108     Int_t nRecPoints = 0;\r
 109     fNRecPoints = 0;\r
 110     fUsedSize = 0;\r
 111     fNDigits = nDigits;\r
 112     fRecPointDataPtr = fFirstRecPointPtr;\r
 113 \r
 114     // Sort our digits\r
 115     SortDigits();\r
 116 \r
 117     //Clusterization starts\r
 118     for (Int_t i = 0; i < nDigits; i++)\r
 119     {\r
 120         fDigitsInCluster = 0;\r
 121 \r
 122          HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);\r
 123         \r
 124         if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)\r
 125         {\r
 126            // Since we have sorted by energy the next digit will have even lower energy, so we return \r
 127            return fNRecPoints;\r
 128         }\r
 129 \r
 130         if(fDigitsPointerArray[i]->fAssociatedCluster != -1)\r
 131         {\r
 132            // The digit is added to a previous cluster, continue\r
 133            continue;\r
 134         }\r
 135 \r
 136         CheckArray();\r
 137         CheckBuffer();\r
 138 \r
 139         // First digit is placed at the fDigits member variable in the recpoint\r
 140         fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
 141 \r
 142         fRecPointDataPtr->fAmp = 0;\r
 143         fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
 144 \r
 145         // Assigning the digit to this rec point\r
 146         fRecPointDataPtr->fDigits = i;\r
 147         fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
 148 \r
 149         // Incrementing the pointer to be ready for new entry\r
 150         fDigitIndexPtr++;\r
 151 \r
 152         fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
 153     \r
 154         \r
 155         //fDigitsPointerArray[i]->fEnergy = 0;\r
 156         fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;\r
 157         \r
 158         \r
 159         fDigitsInCluster++;\r
 160         nRecPoints++;\r
 161 \r
 162         // Scanning for the neighbours\r
 163         if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
 164         {\r
 165             return -1;\r
 166         }\r
 167 \r
 168         //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);\r
 169 \r
 170         fRecPointDataPtr->fMultiplicity = fDigitsInCluster;\r
 171         fRecPointArray[fNRecPoints] = fRecPointDataPtr;\r
 172 \r
 173         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
 174 \r
 175         fNRecPoints++;\r
 176 \r
 177     }//end of clusterization\r
 178 \r
 179     return nRecPoints;\r
 180 }\r
 181 \r
 182 Int_t\r
 183 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
 184 {\r
 185     //see header file for documentation\r
 186 \r
 187     // The following cuts can be used if we sort by posisiton. Not tested, but it should be fine...\r
 188     Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
 189     Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
 190 \r
 191     // All digits for now\r
 192     max = fNDigits;\r
 193     min = 0;\r
 194 \r
 195     for (Int_t j = min; j < max; j++)\r
 196     {\r
 197         if (fDigitsPointerArray[j]->fAssociatedCluster == -1 &&  fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
 198         {\r
 199             if (j != index)\r
 200             {\r
 201                 if (AreNeighbours(fDigitsPointerArray[index],\r
 202                                   fDigitsPointerArray[j]))\r
 203                 {\r
 204                     // Check that the buffer is large enough for adding a digit (can be heavily improved wrt performance)\r
 205                     CheckBuffer();\r
 206 \r
 207                     // Assigning index to digit\r
 208                     *fDigitIndexPtr = j;\r
 209                     fUsedSize += sizeof(Int_t);\r
 210 \r
 211                     // Incrementing digit pointer to be ready for new entry\r
 212                     fDigitIndexPtr++;\r
 213 \r
 214                     // Adding the digit energy to the rec point\r
 215                     fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;\r
 216 \r
 217                     // Setting energy to 0\r
 218                     //fDigitsPointerArray[j]->fEnergy = 0;\r
 219                     \r
 220                     // Setting the associated cluster \r
 221                     fDigitsPointerArray[j]->fAssociatedCluster = fNRecPoints;\r
 222                     \r
 223                     HLTDebug("Added digit with index: %d, energy: %f, to associated cluster: %d", fDigitsPointerArray[j]->fID, fDigitsPointerArray[j]->fEnergy, fDigitsPointerArray[j]->fAssociatedCluster);\r
 224                     \r
 225                     fDigitsInCluster++;\r
 226 \r
 227                     // Scan for neighbours of this digit\r
 228                     ScanForNeighbourDigits(j, recPoint);\r
 229                 }\r
 230             }\r
 231         }\r
 232     }\r
 233     return 0;\r
 234 }\r
 235 \r
 236 Int_t\r
 237 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,\r
 238                                      AliHLTCaloDigitDataStruct* digit2)\r
 239 {\r
 240     //see header file for documentation\r
 241     if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
 242     {\r
 243         Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );\r
 244         Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );\r
 245 \r
 246         // As in the offline code we define neighbours as cells that share an edge, a corner is not  enough\r
 247         //        if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
 248         if (( coldiff <= 1) || ( rowdiff <= 1 ))\r
 249         {\r
 250             // Check also for time\r
 251             if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
 252             {\r
 253                 return 1;\r
 254             }\r
 255         }\r
 256     }\r
 257     return 0;\r
 258 }\r
 259 \r
 260 \r
 261 \r
 262 Int_t AliHLTCaloClusterizer::CheckArray()\r
 263 {\r
 264     // See header file for class documentation\r
 265     if (fArraySize == fNRecPoints)\r
 266     {\r
 267         fArraySize *= 2;\r
 268         AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
 269         memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
 270         delete [] fRecPointArray;\r
 271         fRecPointArray = tmp;\r
 272     }\r
 273     return 0;\r
 274 }\r
 275 \r
 276 Int_t AliHLTCaloClusterizer::CheckBuffer()\r
 277 {\r
 278     // See header file for class documentation\r
 279     if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))\r
 280     {\r
 281         Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
 282         Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);\r
 283         UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
 284         \r
 285         if (tmp == NULL)\r
 286           {\r
 287             HLTError("Pointer error");\r
 288             return(-1);\r
 289           }\r
 290         \r
 291         memcpy(tmp, fFirstRecPointPtr, fUsedSize);\r
 292         fAvailableSize *= 2;\r
 293         for (Int_t n = 0; n < fNRecPoints; n++)\r
 294         {\r
 295             fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
 296         }\r
 297         delete [] fBuffer; //FR\r
 298         fBuffer = tmp; //FR\r
 299         \r
 300         fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
 301         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
 302         fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);\r
 303         //fUsedSize = 0;\r
 304     }\r
 305     return 0;\r
 306 }\r
 307 \r
 308 void AliHLTCaloClusterizer::SetSortDigitsByPosition()\r
 309 {\r
 310     // Sort the digit pointers by position\r
 311     fCompareFunction = &CompareDigitsByPosition;\r
 312     fSortDigits = true;\r
 313     fSortedByPosition = true;\r
 314 }\r
 315 \r
 316 void AliHLTCaloClusterizer::SetSortDigitsByEnergy()\r
 317 {\r
 318     // See header file for class documentation\r
 319     fCompareFunction = &CompareDigitsByEnergy;\r
 320     fSortDigits = true;\r
 321     fSortedByEnergy = true;\r
 322 }\r
 323 \r
 324 void AliHLTCaloClusterizer::SortDigits()\r
 325 {\r
 326     // See header file for class documentation\r
 327     if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);\r
 328 }\r
 329 \r
 330 Int_t\r
 331 AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)\r
 332 {\r
 333     // See header file for documentation\r
 334     return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;\r
 335 }\r
 336 \r
 337 Int_t\r
 338 AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)\r
 339 {\r
 340     // See header file for documentation\r
 341   if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;\r
 342   return 1;\r
 343 }\r