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   fRecPointArray(0),\r
  44   fRecPointDataPtr(0),\r
  45   fFirstRecPointPtr(0),\r
  46   fArraySize(0),\r
  47   fAvailableSize(0),\r
  48   fUsedSize(0),\r
  49   fNRecPoints(0),\r
  50   fDigitIndexPtr(0),\r
  51   fEmcClusteringThreshold(0),\r
  52   fEmcMinEnergyThreshold(0),\r
  53   fEmcTimeGate(0),\r
  54   fDigitsInCluster(0),\r
  55   fDigitsPointerArray(0),\r
  56   fDigitContainerPtr(0),\r
  57   fMaxDigitIndexDiff(0),\r
  58   fNDigits(0)\r
  59 {\r
  60   //See header file for documentation\r
  61   fEmcClusteringThreshold = 0.2;\r
  62   fEmcMinEnergyThreshold = 0.03;\r
  63   fEmcTimeGate = 1.e-6 ;\r
  64   \r
  65   fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
  66   \r
  67   \r
  68   fArraySize = 10;\r
  69   fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
  70   \r
  71   fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
  72   fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(new UChar_t[fAvailableSize]);\r
  73   fFirstRecPointPtr = fRecPointDataPtr;  \r
  74   printf("Start of rec point data: %x, end of rec point data: %x\n", fRecPointDataPtr, reinterpret_cast<UChar_t*>(fRecPointDataPtr) + fAvailableSize);\r
  75 \r
  76 }//end\r
  77 \r
  78 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()  \r
  79 {\r
  80   //See header file for documentation\r
  81 }\r
  82 \r
  83 void \r
  84 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
  85 {\r
  86   // See header file for documentation\r
  87   fRecPointDataPtr = recPointDataPtr;\r
  88 }\r
  89 \r
  90 Int_t \r
  91 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
  92 {\r
  93   //see header file for documentation\r
  94   Int_t nRecPoints = 0;\r
  95   fNRecPoints = 0;\r
  96   fUsedSize = 0;\r
  97   fNDigits = nDigits;\r
  98   fRecPointDataPtr = fFirstRecPointPtr;\r
  99   //Clusterization starts\r
 100   for(Int_t i = 0; i < nDigits; i++)\r
 101     { \r
 102       fDigitsInCluster = 0;\r
 103       //      printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
 104       if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
 105         {\r
 106           continue;\r
 107         }\r
 108       CheckArray();\r
 109       CheckBuffer();\r
 110       //            printf("cluster candidate!\n");\r
 111       // First digit is placed at the fDigits member variable in the recpoint\r
 112       fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
 113 \r
 114       fRecPointDataPtr->fAmp = 0;\r
 115       fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
 116 \r
 117       // Assigning the digit to this rec point\r
 118       fRecPointDataPtr->fDigits = i;\r
 119       printf("Clusterizier: adding digit:  index pointer: %x, index: %d\n", fDigitIndexPtr, *fDigitIndexPtr);\r
 120       fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
 121       \r
 122       // Incrementing the pointer to be ready for new entry\r
 123       fDigitIndexPtr++;\r
 124 \r
 125       fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
 126       fDigitsPointerArray[i]->fEnergy = 0;\r
 127       fDigitsInCluster++;\r
 128       nRecPoints++;\r
 129 \r
 130       // Scanning for the neighbours\r
 131       if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
 132       {\r
 133          return -1;\r
 134       }\r
 135 \r
 136       //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);   \r
 137       \r
 138       fRecPointDataPtr->fMultiplicity = fDigitsInCluster;     \r
 139       printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
 140       printf("Multiplicity: %d\n", fDigitsInCluster);\r
 141       fRecPointArray[fNRecPoints] = fRecPointDataPtr; \r
 142       \r
 143       fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
 144       \r
 145     }//end of clusterization\r
 146    fNRecPoints = nRecPoints;\r
 147    return nRecPoints;\r
 148 }\r
 149 \r
 150 Int_t\r
 151 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
 152 {\r
 153   //see header file for documentation\r
 154   Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
 155   Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
 156 \r
 157   max = fNDigits;\r
 158   min = 0;\r
 159   for(Int_t j = min; j < max; j++)\r
 160     {\r
 161       if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
 162         {\r
 163           if(j != index)\r
 164             {\r
 165               if(AreNeighbours(fDigitsPointerArray[index],\r
 166                                fDigitsPointerArray[j]))\r
 167                 {\r
 168 //                if((fAvailableSize - fUsedSize) < sizeof(Int_t))\r
 169 //                   {\r
 170 //                      UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
 171 //                      memcpy(tmp, fRecPointDataPtr, fAvailableSize);\r
 172 //                      for(Int_t n = 0; n < fNRecPoints; n++)\r
 173 //                      {\r
 174 //                         fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
 175 //                      }\r
 176 //                      fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
 177 //                      fFirstRecPointPtr = fRecPointDataPtr;\r
 178 //                      fUsedSize = 0;\r
 179 //                   }  \r
 180                   CheckBuffer();\r
 181                   // Assigning index to digit\r
 182                   printf("Digit index pointer: %x\n", fDigitIndexPtr);\r
 183                   *fDigitIndexPtr = j;\r
 184                   fUsedSize += sizeof(Int_t);\r
 185                   \r
 186                   printf("Clusterizier: adding digit:  index pointer: %x, index: %d\n", fDigitIndexPtr, *fDigitIndexPtr); \r
 187                   // Incrementing digit pointer to be ready for new entry\r
 188                   fDigitIndexPtr++;\r
 189 \r
 190                   recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
 191                   fDigitsPointerArray[j]->fEnergy = 0;        \r
 192                   fDigitsInCluster++;\r
 193                   ScanForNeighbourDigits(j, recPoint);\r
 194                 }\r
 195             }\r
 196         }\r
 197     }\r
 198   return 0;\r
 199 }\r
 200 \r
 201 Int_t \r
 202 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
 203                                             AliHLTCaloDigitDataStruct* digit2)\r
 204 {\r
 205   //see header file for documentation\r
 206   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
 207     { \r
 208       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
 209       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
 210       if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
 211         {\r
 212 //        cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 213 //          " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 214 \r
 215           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
 216             {\r
 217               return 1; \r
 218             }\r
 219         }\r
 220 \r
 221    /*   Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
 222       Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
 223       if (( coldiff <= 2.4   &&  rowdiff < 0.4 ) || ( coldiff < 0.4 &&  rowdiff <= 2.4 ))\r
 224         {\r
 225           //      cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 226           //        " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 227 \r
 228           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
 229             {\r
 230               return 1; \r
 231             }\r
 232         }\r
 233    */\r
 234       else\r
 235         {\r
 236 //          cout << "Not neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 237 //          " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 238         }\r
 239     }\r
 240   return 0;\r
 241 }\r
 242 \r
 243 \r
 244 \r
 245 Int_t AliHLTCaloClusterizer::CheckArray()\r
 246 {\r
 247       printf("CheckArray: fArraySize: %d, fNRecPoints: %d\n", fArraySize, fNRecPoints);\r
 248       if(fArraySize == fNRecPoints)\r
 249         {\r
 250            printf("Expanding array...");\r
 251            fArraySize *= 2;\r
 252            AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
 253            memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
 254            delete fRecPointArray;\r
 255            fRecPointArray = tmp;\r
 256         }\r
 257    return 0;\r
 258 }\r
 259 \r
 260 Int_t AliHLTCaloClusterizer::CheckBuffer()\r
 261 {\r
 262    // See header file for class documentation \r
 263          printf("CheckBuffer: Used size %d, fAvailableSize: %d\n", fUsedSize, fAvailableSize);\r
 264         if((fAvailableSize - fUsedSize) < sizeof(AliHLTCaloRecPointDataStruct) )\r
 265         {\r
 266            printf("Expanding buffer...\n");\r
 267             Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
 268             fAvailableSize *= 2;\r
 269             UChar_t *tmp = new UChar_t[fAvailableSize];\r
 270             memcpy(tmp, fRecPointDataPtr, fAvailableSize/2);\r
 271             for(Int_t n = 0; n < fNRecPoints; n++)\r
 272             {\r
 273                fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
 274             }\r
 275             delete fRecPointDataPtr;\r
 276             fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
 277             fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
 278             fUsedSize = 0;\r
 279         }\r
 280    return 0;\r
 281 }\r
 282 \r
 283 Int_t AliHLTCaloClusterizer::CheckDigits(AliHLTCaloRecPointDataStruct** recArray, AliHLTCaloDigitDataStruct** digitArray, Int_t nRP)\r
 284 {\r
 285   AliHLTCaloRecPointDataStruct **recpoints = recArray;\r
 286   AliHLTCaloDigitDataStruct **digits = digitArray;\r
 287   Int_t nRecPoints = nRP;\r
 288   \r
 289   if(recArray == 0)\r
 290   {\r
 291      recpoints = fRecPointArray;\r
 292   }\r
 293   if(digitArray == 0)\r
 294   {\r
 295      digits = fDigitsPointerArray;\r
 296   }\r
 297   if(nRP == 0)\r
 298   {\r
 299      nRecPoints = fNRecPoints;\r
 300   }\r
 301   printf("CL: CheckDigits: Number of rec points: %d\n", nRecPoints);\r
 302   for(Int_t i = 0; i < nRecPoints; i++)\r
 303   {\r
 304           \r
 305      AliHLTCaloRecPointDataStruct *recPoint = recpoints[i];\r
 306 \r
 307      //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[0];\r
 308      Int_t multiplicity = recPoint->fMultiplicity;\r
 309      Int_t *digitIndexPtr = &(recPoint->fDigits);\r
 310      printf("CL: Rec point with energy: %f, multiplicity: %d\n", recPoint->fAmp, recPoint->fMultiplicity);\r
 311      for(Int_t j = 0; j < multiplicity; j++)\r
 312      {\r
 313         //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[j];\r
 314         AliHLTCaloDigitDataStruct *digit = digits[*digitIndexPtr];\r
 315         printf("CL: Digit ID: %d, energy: %f, index: %d, indexpointer: %x\n", digit->fID, digit->fEnergy, *digitIndexPtr, digitIndexPtr);\r
 316         digitIndexPtr++;\r
 317         //recPoint = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitIndexPtr);\r
 318      }\r
 319   }\r
 320      \r
 321      \r
 322      \r
 323      \r
 324      \r
 325      \r
 326      \r
 327 }\r
 328 \r
 329 Int_t AliHLTCaloClusterizer::CheckDigits(AliHLTCaloRecPointDataStruct** recArray, AliHLTCaloDigitDataStruct* digitArray, Int_t nRP)\r
 330 {\r
 331   AliHLTCaloRecPointDataStruct **recpoints = recArray;\r
 332   AliHLTCaloDigitDataStruct *digits = digitArray;\r
 333   Int_t nRecPoints = nRP;\r
 334   \r
 335   if(recArray == 0)\r
 336   {\r
 337      recpoints = fRecPointArray;\r
 338   }\r
 339     if(nRP == 0)\r
 340   {\r
 341      nRecPoints = fNRecPoints;\r
 342   }\r
 343   printf("CL: CheckDigits: Number of rec points: %d\n", nRecPoints);\r
 344   for(Int_t i = 0; i < nRecPoints; i++)\r
 345   {\r
 346           \r
 347      AliHLTCaloRecPointDataStruct *recPoint = recpoints[i];\r
 348 \r
 349      //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[0];\r
 350      Int_t multiplicity = recPoint->fMultiplicity;\r
 351      Int_t *digitIndexPtr = &(recPoint->fDigits);\r
 352      printf("CL: Rec point with energy: %f, multiplicity: %d\n", recPoint->fAmp, recPoint->fMultiplicity);\r
 353      for(Int_t j = 0; j < multiplicity; j++)\r
 354      {\r
 355         //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[j];\r
 356         AliHLTCaloDigitDataStruct digit = digits[*digitIndexPtr];\r
 357         printf("CL: digits: %x, recpoints: %x, digitIndexPtr: %x\n", digits, recpoints, digitIndexPtr);\r
 358         printf("CL: Digit ID: %d, energy: %f, index: %d, indexpointer: %x\n", digit.fID, digit.fEnergy, *digitIndexPtr, digitIndexPtr);\r
 359         digitIndexPtr++;\r
 360         //recPoint = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitIndexPtr);\r
 361      }\r
 362   }\r
 363      \r
 364      \r
 365      \r
 366      \r
 367      \r
 368      \r
 369      \r
 370 }