[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterizer.cxx

// $Id$\r
\r
/**************************************************************************\r
 * This file is property of and copyright by the ALICE HLT Project        * \r
 * All rights reserved.                                                   *\r
 *                                                                        *\r
 * Primary Authors: Oystein Djuvsland                                     *\r
 *                                                                        *\r
 * Permission to use, copy, modify and distribute this software and its   *\r
 * documentation strictly for non-commercial purposes is hereby granted   *\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
 * provided "as is" without express or implied warranty.                  *\r
 **************************************************************************/\r
\r
/** \r
 * @file   AliHLTCaloClusterizer.cxx\r
 * @author Oystein Djuvsland\r
 * @date \r
 * @brief  Clusterizer for PHOS HLT \r
 */\r
\r
// see header file for class documentation\r
// or\r
// refer to README to build package\r
// or\r
// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt\r
\r
#include "AliHLTCaloClusterizer.h"\r
#include "AliHLTLogging.h"\r
#include "TMath.h"\r
#include "AliHLTCaloRecPointDataStruct.h"\r
#include "AliHLTCaloDigitDataStruct.h"\r
#include "AliHLTCaloDigitContainerDataStruct.h"\r
#include "AliHLTCaloConstantsHandler.h"\r
\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
{\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
\r
}//end\r
\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
\r
AliHLTCaloClusterizer::~AliHLTCaloClusterizer()  \r
{\r
  //See header file for documentation\r
}\r
\r
void \r
AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
{\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
{\r
  //see header file for documentation\r
  Int_t nRecPoints = 0;\r
  \r
  fAvailableSize = availableSize;\r
  \r
  fNDigits = nDigits;\r
\r
  UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \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
	\r
	if(fAvailableSize < (sizeof(AliHLTCaloRecPointDataStruct)))\r
	{\r
	  HLTError("Out of buffer, stopping clusterisation");\r
	  return -1; \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
\r
      // Assigning the digit to this rec point\r
      fRecPointDataPtr->fDigits = i;\r
\r
      // Incrementing the pointer to be ready for new entry\r
      fDigitIndexPtr++;\r
\r
      fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
      fDigitsPointerArray[i]->fEnergy = 0;\r
      fDigitsInCluster++;\r
      nRecPoints++;\r
\r
      // Scanning for the neighbours\r
      if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
      {\r
	 return -1;\r
      }\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
\r
    }//end of clusterization\r
\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
\r
  max = fNDigits;\r
  min = 0;\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
    }\r
  return 0;\r
}\r
\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
\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
    }\r
  return 0;\r
}\r
Commit	Line	Data
ef44ec64	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
ef44ec64	32	#include "AliHLTLogging.h"\r
ef44ec64	33	#include "TMath.h"\r
ef44ec64	34	#include "AliHLTCaloRecPointDataStruct.h"\r
	35	#include "AliHLTCaloDigitDataStruct.h"\r
	36	#include "AliHLTCaloDigitContainerDataStruct.h"\r
4f4b7ba4	37	#include "AliHLTCaloConstantsHandler.h"\r
ef44ec64	38	\r
	39	ClassImp(AliHLTCaloClusterizer);\r
	40	\r
	41	AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):\r
4f4b7ba4	42	AliHLTCaloConstantsHandler(det),\r
ef44ec64	43	fRecPointDataPtr(0),\r
ad44d760	44	fDigitIndexPtr(0),\r
ef44ec64	45	fEmcClusteringThreshold(0),\r
	46	fEmcMinEnergyThreshold(0),\r
	47	fEmcTimeGate(0),\r
	48	fDigitsInCluster(0),\r
f92dcf21	49	fDigitsPointerArray(0),\r
ef44ec64	50	fDigitContainerPtr(0),\r
ad44d760	51	fMaxDigitIndexDiff(0),\r
ad44d760	52	fNDigits(0)\r
ef44ec64	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
c375e15d	60	\r
c375e15d	61	}//end\r
c375e15d	62	\r
4f4b7ba4	63	\r
4f4b7ba4	64	//BALLE how do you set the right detector?\r
ad44d760	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
ef44ec64	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
ef44ec64	92	Int_t \r
ad44d760	93	AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits, UInt_t availableSize, UInt_t& totSize)\r
ef44ec64	94	{\r
	95	//see header file for documentation\r
	96	Int_t nRecPoints = 0;\r
9bd6f29b	97	\r
	98	fAvailableSize = availableSize;\r
	99	\r
ad44d760	100	fNDigits = nDigits;\r
ad44d760	101	\r
ef44ec64	102	UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \r
	103	\r
	104	//Clusterization starts\r
f92dcf21	105	for(Int_t i = 0; i < nDigits; i++)\r
ef44ec64	106	{ \r
ef44ec64	107	fDigitsInCluster = 0;\r
ad44d760	108	// printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
ad44d760	109	if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
ef44ec64	110	{\r
	111	continue;\r
	112	}\r
9bd6f29b	113	\r
9bd6f29b	114	if(fAvailableSize < (sizeof(AliHLTCaloRecPointDataStruct)))\r
ef44ec64	115	{\r
9bd6f29b	116	HLTError("Out of buffer, stopping clusterisation");\r
9bd6f29b	117	return -1; \r
ef44ec64	118	}\r
9bd6f29b	119	\r
ad44d760	120	// printf("cluster candidate!\n");\r
ef44ec64	121	// First digit is placed at the fDigits member variable in the recpoint\r
ad44d760	122	fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
ef44ec64	123	\r
ef44ec64	124	fRecPointDataPtr->fAmp = 0;\r
ad44d760	125	fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
ef44ec64	126	\r
ad44d760	127	// Assigning the digit to this rec point\r
ad44d760	128	fRecPointDataPtr->fDigits = i;\r
ef44ec64	129	\r
ef44ec64	130	// Incrementing the pointer to be ready for new entry\r
ad44d760	131	fDigitIndexPtr++;\r
ef44ec64	132	\r
ad44d760	133	fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
ad44d760	134	fDigitsPointerArray[i]->fEnergy = 0;\r
ef44ec64	135	fDigitsInCluster++;\r
	136	nRecPoints++;\r
	137	\r
	138	// Scanning for the neighbours\r
9bd6f29b	139	if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
	140	{\r
	141	return -1;\r
	142	}\r
ef44ec64	143	\r
	144	totSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct); \r
	145	fRecPointDataPtr->fMultiplicity = fDigitsInCluster; \r
ad44d760	146	// printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
ad44d760	147	fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
ef44ec64	148	\r
ef44ec64	149	}//end of clusterization\r
	150	\r
	151	return nRecPoints;\r
	152	}\r
	153	\r
9bd6f29b	154	Int_t\r
ef44ec64	155	AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
	156	{\r
	157	//see header file for documentation\r
ad44d760	158	Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
ef44ec64	159	Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
ef44ec64	160	\r
ad44d760	161	max = fNDigits;\r
ef44ec64	162	min = 0;\r
	163	for(Int_t j = min; j < max; j++)\r
	164	{\r
ad44d760	165	if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
ef44ec64	166	{\r
	167	if(j != index)\r
	168	{\r
ad44d760	169	if(AreNeighbours(fDigitsPointerArray[index],\r
ad44d760	170	fDigitsPointerArray[j]))\r
ef44ec64	171	{\r
9bd6f29b	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
ad44d760	178	*fDigitIndexPtr = j;\r
ef44ec64	179	// Incrementing digit pointer to be ready for new entry\r
ad44d760	180	fDigitIndexPtr++;\r
ef44ec64	181	\r
ad44d760	182	recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
ad44d760	183	fDigitsPointerArray[j]->fEnergy = 0; \r
ef44ec64	184	fDigitsInCluster++;\r
	185	ScanForNeighbourDigits(j, recPoint);\r
	186	}\r
	187	}\r
	188	}\r
	189	}\r
9bd6f29b	190	return 0;\r
ef44ec64	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