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