[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
  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
      if(availableSize < (totSize + maxRecPointSize)) \r
	{\r
	  return -1; //Might get out of buffer, exiting\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
      ScanForNeighbourDigits(i, 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
\r
    }//end of clusterization\r
\r
   return nRecPoints;\r
}\r
\r
void\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
		  // Assigning value to digit ptr\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;\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
	97	\r
ad44d760	98	fNDigits = nDigits;\r
ad44d760	99	\r
ef44ec64	100	UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \r
	101	\r
	102	//Clusterization starts\r
f92dcf21	103	for(Int_t i = 0; i < nDigits; i++)\r
ef44ec64	104	{ \r
ef44ec64	105	fDigitsInCluster = 0;\r
ad44d760	106	// printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
ad44d760	107	if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
ef44ec64	108	{\r
	109	continue;\r
	110	}\r
	111	if(availableSize < (totSize + maxRecPointSize)) \r
	112	{\r
	113	return -1; //Might get out of buffer, exiting\r
	114	}\r
ad44d760	115	// printf("cluster candidate!\n");\r
ef44ec64	116	// First digit is placed at the fDigits member variable in the recpoint\r
ad44d760	117	fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
ef44ec64	118	\r
ef44ec64	119	fRecPointDataPtr->fAmp = 0;\r
ad44d760	120	fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
ef44ec64	121	\r
ad44d760	122	// Assigning the digit to this rec point\r
ad44d760	123	fRecPointDataPtr->fDigits = i;\r
ef44ec64	124	\r
ef44ec64	125	// Incrementing the pointer to be ready for new entry\r
ad44d760	126	fDigitIndexPtr++;\r
ef44ec64	127	\r
ad44d760	128	fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
ad44d760	129	fDigitsPointerArray[i]->fEnergy = 0;\r
ef44ec64	130	fDigitsInCluster++;\r
	131	nRecPoints++;\r
	132	\r
	133	// Scanning for the neighbours\r
	134	ScanForNeighbourDigits(i, fRecPointDataPtr);\r
	135	\r
	136	totSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct); \r
	137	fRecPointDataPtr->fMultiplicity = fDigitsInCluster; \r
ad44d760	138	// printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
ad44d760	139	fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
ef44ec64	140	\r
ef44ec64	141	}//end of clusterization\r
	142	\r
	143	return nRecPoints;\r
	144	}\r
	145	\r
	146	void\r
	147	AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
	148	{\r
	149	//see header file for documentation\r
ad44d760	150	Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
ef44ec64	151	Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
ef44ec64	152	\r
ad44d760	153	max = fNDigits;\r
ef44ec64	154	min = 0;\r
	155	for(Int_t j = min; j < max; j++)\r
	156	{\r
ad44d760	157	if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
ef44ec64	158	{\r
	159	if(j != index)\r
	160	{\r
ad44d760	161	if(AreNeighbours(fDigitsPointerArray[index],\r
ad44d760	162	fDigitsPointerArray[j]))\r
ef44ec64	163	{\r
ef44ec64	164	// Assigning value to digit ptr\r
ad44d760	165	*fDigitIndexPtr = j;\r
ef44ec64	166	// Incrementing digit pointer to be ready for new entry\r
ad44d760	167	fDigitIndexPtr++;\r
ef44ec64	168	\r
ad44d760	169	recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
ad44d760	170	fDigitsPointerArray[j]->fEnergy = 0; \r
ef44ec64	171	fDigitsInCluster++;\r
	172	ScanForNeighbourDigits(j, recPoint);\r
	173	}\r
	174	}\r
	175	}\r
	176	}\r
	177	return;\r
	178	}\r
	179	\r
	180	Int_t \r
	181	AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
	182	AliHLTCaloDigitDataStruct* digit2)\r
	183	{\r
	184	//see header file for documentation\r
	185	if ( (digit1->fModule == digit2->fModule) /&& (coord1[1]==coord2[1])/ ) // inside the same PHOS module\r
	186	{ \r
	187	// Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ ); \r
	188	// Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
	189	// if (( coldiff <= 1 && rowdiff == 0 ) \|\| ( coldiff == 0 && rowdiff <= 1 ))\r
	190	// {\r
	191	// cout << "Are neighbours: digit (E = " << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
	192	// " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
	193	\r
	194	// if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
	195	// {\r
	196	// return 1; \r
	197	// }\r
	198	// }\r
	199	\r
	200	Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ ); \r
	201	Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
	202	if (( coldiff <= 2.4 && rowdiff < 0.4 ) \|\| ( coldiff < 0.4 && rowdiff <= 2.4 ))\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	else\r
	213	{\r
	214	// cout << "Not neighbours: digit (E = " << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
	215	// " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
	216	}\r
	217	}\r
	218	return 0;\r
	219	}\r