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   fRecPointDataPtr(0),\r
  44   fDigitIndexPtr(0),\r
  45   fEmcClusteringThreshold(0),\r
  46   fEmcMinEnergyThreshold(0),\r
  47   fEmcTimeGate(0),\r
  48   fDigitsInCluster(0),\r
  49   fDigitsPointerArray(0),\r
  50   fDigitContainerPtr(0),\r
  51   fMaxDigitIndexDiff(0),\r
  52   fNDigits(0)\r
  53 {\r
  54   //See header file for documentation\r
  55   fEmcClusteringThreshold = 0.2;\r
  56   fEmcMinEnergyThreshold = 0.03;\r
  57   fEmcTimeGate = 1.e-6 ;\r
  58   \r
  59   fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
  60 \r
  61 }//end\r
  62 \r
  63 \r
  64 //BALLE how do you set the right detector?\r
  65 // AliHLTCaloClusterizer::AliHLTCaloClusterizer(const AliHLTCaloClusterizer &) :\r
  66 //   AliHLTCaloConstantsHandler("BALLE"),\r
  67 //   fRecPointDataPtr(0),\r
  68 //   fDigitDataPtr(0),\r
  69 //   fEmcClusteringThreshold(0),\r
  70 //   fEmcMinEnergyThreshold(0),\r
  71 //   fEmcTimeGate(0),\r
  72 //   fDigitsInCluster(0),\r
  73 //   fDigitContainerPtr(0),\r
  74 //   fMaxDigitIndexDiff(0)\r
  75 // {\r
  76 //   // dummy copy constructor\r
  77 // }//end\r
  78 \r
  79 \r
  80 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()  \r
  81 {\r
  82   //See header file for documentation\r
  83 }\r
  84 \r
  85 void \r
  86 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
  87 {\r
  88   // See header file for documentation\r
  89   fRecPointDataPtr = recPointDataPtr;\r
  90 }\r
  91 \r
  92 Int_t \r
  93 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits, UInt_t availableSize, UInt_t& totSize)\r
  94 {\r
  95   //see header file for documentation\r
  96   Int_t nRecPoints = 0;\r
  97   \r
  98   fAvailableSize = availableSize;\r
  99   \r
 100   fNDigits = nDigits;\r
 101 \r
 102   UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \r
 103 \r
 104   //Clusterization starts\r
 105   for(Int_t i = 0; i < nDigits; i++)\r
 106     { \r
 107       fDigitsInCluster = 0;\r
 108       //      printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
 109       if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
 110         {\r
 111           continue;\r
 112         }\r
 113         \r
 114         if(fAvailableSize < (sizeof(AliHLTCaloRecPointDataStruct)))\r
 115         {\r
 116           HLTError("Out of buffer, stopping clusterisation");\r
 117           return -1; \r
 118         }\r
 119         \r
 120       //            printf("cluster candidate!\n");\r
 121       // First digit is placed at the fDigits member variable in the recpoint\r
 122       fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
 123 \r
 124       fRecPointDataPtr->fAmp = 0;\r
 125       fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
 126 \r
 127       // Assigning the digit to this rec point\r
 128       fRecPointDataPtr->fDigits = i;\r
 129 \r
 130       // Incrementing the pointer to be ready for new entry\r
 131       fDigitIndexPtr++;\r
 132 \r
 133       fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
 134       fDigitsPointerArray[i]->fEnergy = 0;\r
 135       fDigitsInCluster++;\r
 136       nRecPoints++;\r
 137 \r
 138       // Scanning for the neighbours\r
 139       if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
 140       {\r
 141          return -1;\r
 142       }\r
 143 \r
 144       totSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);   \r
 145       fRecPointDataPtr->fMultiplicity = fDigitsInCluster;     \r
 146       //      printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
 147       fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
 148 \r
 149     }//end of clusterization\r
 150 \r
 151    return nRecPoints;\r
 152 }\r
 153 \r
 154 Int_t\r
 155 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
 156 {\r
 157   //see header file for documentation\r
 158   Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
 159   Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
 160 \r
 161   max = fNDigits;\r
 162   min = 0;\r
 163   for(Int_t j = min; j < max; j++)\r
 164     {\r
 165       if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
 166         {\r
 167           if(j != index)\r
 168             {\r
 169               if(AreNeighbours(fDigitsPointerArray[index],\r
 170                                fDigitsPointerArray[j]))\r
 171                 {\r
 172                   if(fAvailableSize < (sizeof(Int_t)))\r
 173                      {\r
 174                         HLTError("Out of buffer, stopping clusterisation");\r
 175                         return -1; \r
 176                      }  \r
 177                   // Assigning index to digit\r
 178                   *fDigitIndexPtr = j;\r
 179                   // Incrementing digit pointer to be ready for new entry\r
 180                   fDigitIndexPtr++;\r
 181 \r
 182                   recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
 183                   fDigitsPointerArray[j]->fEnergy = 0;        \r
 184                   fDigitsInCluster++;\r
 185                   ScanForNeighbourDigits(j, recPoint);\r
 186                 }\r
 187             }\r
 188         }\r
 189     }\r
 190   return 0;\r
 191 }\r
 192 \r
 193 Int_t \r
 194 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
 195                                             AliHLTCaloDigitDataStruct* digit2)\r
 196 {\r
 197   //see header file for documentation\r
 198   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
 199     { \r
 200 //       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
 201 //       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
 202 //       if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
 203 //      {\r
 204 //        cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 205 //          " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 206 \r
 207 //        if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
 208 //          {\r
 209 //            return 1; \r
 210 //          }\r
 211 //      }\r
 212 \r
 213       Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
 214       Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
 215       if (( coldiff <= 2.4   &&  rowdiff < 0.4 ) || ( coldiff < 0.4 &&  rowdiff <= 2.4 ))\r
 216         {\r
 217           //      cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 218           //        " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 219 \r
 220           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
 221             {\r
 222               return 1; \r
 223             }\r
 224         }\r
 225       else\r
 226         {\r
 227           //  cout << "Not neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
 228           //  " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
 229         }\r
 230     }\r
 231   return 0;\r
 232 }\r