// $Id$\r
\r
/**************************************************************************\r
- * This file is property of and copyright by the ALICE HLT Project * \r
+ * This file is property of and copyright by the ALICE HLT Project *\r
* All rights reserved. *\r
* *\r
* Primary Authors: Oystein Djuvsland *\r
* without fee, provided that the above copyright notice appears in all *\r
* copies and that both the copyright notice and this permission notice *\r
* appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is * \r
+ * about the suitability of this software for any purpose. It is *\r
* provided "as is" without express or implied warranty. *\r
**************************************************************************/\r
\r
-/** \r
+/**\r
* @file AliHLTCaloClusterizer.cxx\r
* @author Oystein Djuvsland\r
- * @date \r
- * @brief Clusterizer for PHOS HLT \r
+ * @date\r
+ * @brief Clusterizer for PHOS HLT\r
*/\r
\r
// see header file for class documentation\r
ClassImp(AliHLTCaloClusterizer);\r
\r
AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):\r
- AliHLTCaloConstantsHandler(det),\r
- fRecPointDataPtr(0),\r
- fDigitIndexPtr(0),\r
- fEmcClusteringThreshold(0),\r
- fEmcMinEnergyThreshold(0),\r
- fEmcTimeGate(0),\r
- fDigitsInCluster(0),\r
- fDigitsPointerArray(0),\r
- fDigitContainerPtr(0),\r
- fMaxDigitIndexDiff(0),\r
- fNDigits(0)\r
+ AliHLTCaloConstantsHandler(det),\r
+ fCompareFunction(CompareDigitsByPosition),\r
+ fRecPointArray(0),\r
+ fRecPointDataPtr(0),\r
+ fFirstRecPointPtr(0),\r
+ fArraySize(0),\r
+ fAvailableSize(0),\r
+ fUsedSize(0),\r
+ fNRecPoints(0),\r
+ fDigitIndexPtr(0),\r
+ fEmcClusteringThreshold(0),\r
+ fEmcMinEnergyThreshold(0),\r
+ fEmcTimeGate(0),\r
+ fDigitsInCluster(0),\r
+ fDigitsPointerArray(0),\r
+ fDigitContainerPtr(0),\r
+ fMaxDigitIndexDiff(0),\r
+ fNDigits(0),\r
+ fSortedByPosition(false),\r
+ fSortedByEnergy(false),\r
+ fSortDigits(false),\r
+ fIsEMCAL(false),\r
+ fBuffer(0)\r
+ \r
{\r
- //See header file for documentation\r
- fEmcClusteringThreshold = 0.2;\r
- fEmcMinEnergyThreshold = 0.03;\r
- fEmcTimeGate = 1.e-6 ;\r
- \r
- fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
+ //See header file for documentation\r
+ //fEmcClusteringThreshold = 0.2;\r
+ //fEmcMinEnergyThreshold = 0.03;\r
\r
-}//end\r
+ fEmcClusteringThreshold = 0.1;\r
+ fEmcMinEnergyThreshold = 0.01;\r
+ fEmcTimeGate = 1.e-6 ;\r
\r
+ fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
\r
-//BALLE how do you set the right detector?\r
-// AliHLTCaloClusterizer::AliHLTCaloClusterizer(const AliHLTCaloClusterizer &) :\r
-// AliHLTCaloConstantsHandler("BALLE"),\r
-// fRecPointDataPtr(0),\r
-// fDigitDataPtr(0),\r
-// fEmcClusteringThreshold(0),\r
-// fEmcMinEnergyThreshold(0),\r
-// fEmcTimeGate(0),\r
-// fDigitsInCluster(0),\r
-// fDigitContainerPtr(0),\r
-// fMaxDigitIndexDiff(0)\r
-// {\r
-// // dummy copy constructor\r
-// }//end\r
\r
+ fArraySize = 10;\r
+ fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
+\r
+ fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
+ fBuffer = new UChar_t[fAvailableSize]; //FR\r
+ \r
+ fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);\r
+ fRecPointDataPtr = fFirstRecPointPtr;\r
+\r
+}//end\r
\r
-AliHLTCaloClusterizer::~AliHLTCaloClusterizer() \r
+AliHLTCaloClusterizer::~AliHLTCaloClusterizer()\r
{\r
- //See header file for documentation\r
+ //See header file for documentation\r
+ delete [] fBuffer; //FR\r
+ fBuffer = NULL; //FR\r
+ fAvailableSize = 0; //FR\r
+\r
+ delete [] fRecPointArray;\r
}\r
\r
-void \r
+void\r
AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
{\r
- // See header file for documentation\r
- fRecPointDataPtr = recPointDataPtr;\r
+ // See header file for documentation\r
+ fRecPointDataPtr = recPointDataPtr;\r
}\r
\r
-Int_t \r
-AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits, UInt_t availableSize, UInt_t& totSize)\r
+Int_t\r
+AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
{\r
- //see header file for documentation\r
- Int_t nRecPoints = 0;\r
- \r
- fAvailableSize = availableSize;\r
- \r
- fNDigits = nDigits;\r
+ //see header file for documentation\r
+ Int_t nRecPoints = 0;\r
+ fNRecPoints = 0;\r
+ fUsedSize = 0;\r
+ fNDigits = nDigits;\r
+ fRecPointDataPtr = fFirstRecPointPtr;\r
\r
- UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \r
+ // Sort our digits\r
+ SortDigits();\r
\r
- //Clusterization starts\r
- for(Int_t i = 0; i < nDigits; i++)\r
- { \r
- fDigitsInCluster = 0;\r
- // printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
- if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
- {\r
- continue;\r
- }\r
+ //Clusterization starts\r
+ for (Int_t i = 0; i < nDigits; i++)\r
+ {\r
+ fDigitsInCluster = 0;\r
+\r
+ HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);\r
\r
- if(fAvailableSize < (sizeof(AliHLTCaloRecPointDataStruct)))\r
+ if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)\r
+ {\r
+ // Since we have sorted by energy the next digit will have even lower energy, so we return \r
+ return fNRecPoints;\r
+ }\r
+\r
+ if(fDigitsPointerArray[i]->fAssociatedCluster != -1)\r
{\r
- HLTError("Out of buffer, stopping clusterisation");\r
- return -1; \r
+ // The digit is added to a previous cluster, continue\r
+ continue;\r
}\r
- \r
- // printf("cluster candidate!\n");\r
- // First digit is placed at the fDigits member variable in the recpoint\r
- fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
\r
- fRecPointDataPtr->fAmp = 0;\r
- fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
+ CheckArray();\r
+ CheckBuffer();\r
+\r
+ // First digit is placed at the fDigits member variable in the recpoint\r
+ fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
+\r
+ fRecPointDataPtr->fAmp = 0;\r
+ fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
+\r
+ // Assigning the digit to this rec point\r
+ fRecPointDataPtr->fDigits = i;\r
+ fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
+\r
+ // Incrementing the pointer to be ready for new entry\r
+ fDigitIndexPtr++;\r
\r
- // Assigning the digit to this rec point\r
- fRecPointDataPtr->fDigits = i;\r
+ fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
+ \r
+ \r
+ //fDigitsPointerArray[i]->fEnergy = 0;\r
+ fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;\r
+ \r
+ \r
+ fDigitsInCluster++;\r
+ nRecPoints++;\r
\r
- // Incrementing the pointer to be ready for new entry\r
- fDigitIndexPtr++;\r
+ // Scanning for the neighbours\r
+ if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
+ {\r
+ return -1;\r
+ }\r
\r
- fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
- fDigitsPointerArray[i]->fEnergy = 0;\r
- fDigitsInCluster++;\r
- nRecPoints++;\r
+ //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);\r
\r
- // Scanning for the neighbours\r
- if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
- {\r
- return -1;\r
- }\r
+ fRecPointDataPtr->fMultiplicity = fDigitsInCluster;\r
+ fRecPointArray[fNRecPoints] = fRecPointDataPtr;\r
\r
- totSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct); \r
- fRecPointDataPtr->fMultiplicity = fDigitsInCluster; \r
- // printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
- fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
+ fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
+\r
+ fNRecPoints++;\r
\r
}//end of clusterization\r
\r
- return nRecPoints;\r
+ return nRecPoints;\r
}\r
\r
Int_t\r
AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
{\r
- //see header file for documentation\r
- Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
- Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
+ //see header file for documentation\r
+\r
+ // The following cuts can be used if we sort by posisiton. Not tested, but it should be fine...\r
+ Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
+ Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
+\r
+ // All digits for now\r
+ max = fNDigits;\r
+ min = 0;\r
\r
- max = fNDigits;\r
- min = 0;\r
- for(Int_t j = min; j < max; j++)\r
+ for (Int_t j = min; j < max; j++)\r
{\r
- if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
- {\r
- if(j != index)\r
- {\r
- if(AreNeighbours(fDigitsPointerArray[index],\r
- fDigitsPointerArray[j]))\r
- {\r
- if(fAvailableSize < (sizeof(Int_t)))\r
- {\r
- HLTError("Out of buffer, stopping clusterisation");\r
- return -1; \r
- } \r
- // Assigning index to digit\r
- *fDigitIndexPtr = j;\r
- // Incrementing digit pointer to be ready for new entry\r
- fDigitIndexPtr++;\r
-\r
- recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
- fDigitsPointerArray[j]->fEnergy = 0; \r
- fDigitsInCluster++;\r
- ScanForNeighbourDigits(j, recPoint);\r
- }\r
- }\r
- }\r
+ if (fDigitsPointerArray[j]->fAssociatedCluster == -1 && fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
+ {\r
+ if (j != index)\r
+ {\r
+ if (AreNeighbours(fDigitsPointerArray[index],\r
+ fDigitsPointerArray[j]))\r
+ {\r
+ // Check that the buffer is large enough for adding a digit (can be heavily improved wrt performance)\r
+ CheckBuffer();\r
+\r
+ // Assigning index to digit\r
+ *fDigitIndexPtr = j;\r
+ fUsedSize += sizeof(Int_t);\r
+\r
+ // Incrementing digit pointer to be ready for new entry\r
+ fDigitIndexPtr++;\r
+\r
+ // Adding the digit energy to the rec point\r
+ fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;\r
+\r
+ // Setting energy to 0\r
+ //fDigitsPointerArray[j]->fEnergy = 0;\r
+ \r
+ // Setting the associated cluster \r
+ fDigitsPointerArray[j]->fAssociatedCluster = fNRecPoints;\r
+ \r
+ HLTDebug("Added digit with index: %d, energy: %f, to associated cluster: %d", fDigitsPointerArray[j]->fID, fDigitsPointerArray[j]->fEnergy, fDigitsPointerArray[j]->fAssociatedCluster);\r
+ \r
+ fDigitsInCluster++;\r
+\r
+ // Scan for neighbours of this digit\r
+ ScanForNeighbourDigits(j, recPoint);\r
+ }\r
+ }\r
+ }\r
}\r
- return 0;\r
+ return 0;\r
}\r
\r
-Int_t \r
-AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
- AliHLTCaloDigitDataStruct* digit2)\r
+Int_t\r
+AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,\r
+ AliHLTCaloDigitDataStruct* digit2)\r
{\r
- //see header file for documentation\r
- if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
- { \r
-// Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ ); \r
-// Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
-// if (( coldiff <= 1 && rowdiff == 0 ) || ( coldiff == 0 && rowdiff <= 1 ))\r
-// {\r
-// cout << "Are neighbours: digit (E = " << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
-// " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
-\r
-// if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
-// {\r
-// return 1; \r
-// }\r
-// }\r
-\r
- Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ ); \r
- Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
- if (( coldiff <= 2.4 && rowdiff < 0.4 ) || ( coldiff < 0.4 && rowdiff <= 2.4 ))\r
- {\r
- // cout << "Are neighbours: digit (E = " << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
- // " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
+ //see header file for documentation\r
+ if ( (digit1->fModule == digit2->fModule) || AreEdgeCells(digit1, digit2))\r
+ {\r
+ Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );\r
+ Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );\r
\r
- if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
- {\r
- return 1; \r
- }\r
- }\r
- else\r
- {\r
- // cout << "Not neighbours: digit (E = " << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
- // " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
- }\r
+ // Common edge defines neighbour\r
+ //if (( coldiff <= 1 && rowdiff == 0 ) || ( coldiff == 0 && rowdiff <= 1 ))\r
+ // Common edge and corner defines neighbour\r
+ if (( coldiff <= 1 && rowdiff <= 1 ))\r
+ {\r
+ // Check also for time\r
+ if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
+ {\r
+ return 1;\r
+ }\r
+ }\r
+ }\r
+ return 0;\r
+}\r
+\r
+\r
+\r
+Int_t AliHLTCaloClusterizer::CheckArray()\r
+{\r
+ // See header file for class documentation\r
+ if (fArraySize == fNRecPoints)\r
+ {\r
+ fArraySize *= 2;\r
+ AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
+ memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
+ delete [] fRecPointArray;\r
+ fRecPointArray = tmp;\r
+ }\r
+ return 0;\r
+}\r
+\r
+Int_t AliHLTCaloClusterizer::CheckBuffer()\r
+{\r
+ // See header file for class documentation\r
+ if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))\r
+ {\r
+ Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
+ Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);\r
+ UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
+ \r
+ if (tmp == NULL)\r
+ {\r
+ HLTError("Pointer error");\r
+ return(-1);\r
+ }\r
+ \r
+ memcpy(tmp, fFirstRecPointPtr, fUsedSize);\r
+ fAvailableSize *= 2;\r
+ for (Int_t n = 0; n < fNRecPoints; n++)\r
+ {\r
+ fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
+ }\r
+ delete [] fBuffer; //FR\r
+ fBuffer = tmp; //FR\r
+ \r
+ fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
+ fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
+ fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);\r
+ //fUsedSize = 0;\r
+ }\r
+ return 0;\r
+}\r
+\r
+void AliHLTCaloClusterizer::SetSortDigitsByPosition()\r
+{\r
+ // Sort the digit pointers by position\r
+ fCompareFunction = &CompareDigitsByPosition;\r
+ fSortDigits = true;\r
+ fSortedByPosition = true;\r
+}\r
+\r
+void AliHLTCaloClusterizer::SetSortDigitsByEnergy()\r
+{\r
+ // See header file for class documentation\r
+ fCompareFunction = &CompareDigitsByEnergy;\r
+ fSortDigits = true;\r
+ fSortedByEnergy = true;\r
+}\r
+\r
+void AliHLTCaloClusterizer::SortDigits()\r
+{\r
+ // See header file for class documentation\r
+ if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);\r
+}\r
+\r
+Int_t\r
+AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)\r
+{\r
+ // See header file for documentation\r
+ return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;\r
+}\r
+\r
+Int_t\r
+AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)\r
+{\r
+ // See header file for documentation\r
+ if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;\r
+ return 1;\r
+}\r
+\r
+void AliHLTCaloClusterizer::SetDetector(TString det)\r
+{\r
+ if(det.CompareTo("EMCAL"))\r
+ {\r
+ fIsEMCAL = true;\r
+ }\r
+ else\r
+ {\r
+ fIsEMCAL = false;\r
+ }\r
+}\r
+\r
+Bool_t AliHLTCaloClusterizer::AreEdgeCells(AliHLTCaloDigitDataStruct *digit0, AliHLTCaloDigitDataStruct *digit1)\r
+{\r
+ if(fIsEMCAL)\r
+ {\r
+ Int_t modDiff = digit0->fModule - digit1->fModule;\r
+ if(TMath::Abs(modDiff) > 1) return kFALSE;\r
+ if(digit0->fModule > digit1->fModule && digit1->fModule%2 == 0) \r
+ {\r
+ if(digit0->fZ == 0 && digit1->fZ == (fCaloConstants->GetNZROWSMOD()-1))\r
+ return kTRUE;\r
+ }\r
+ if(digit1->fModule > digit0->fModule && digit0->fModule%2 == 0) \r
+ {\r
+ if(digit1->fZ == 0 && digit0->fZ == (fCaloConstants->GetNZROWSMOD()-1))\r
+ return kTRUE;\r
}\r
- return 0;\r
+ }\r
+ \r
+ return false;\r
+\r
}\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff