HLT/PHOS/AliHLTPHOSClusterizer.cxx

   1 // $Id$
   2
   3 /**************************************************************************
   4  * This file is property of and copyright by the ALICE HLT Project        *
   5  * All rights reserved.                                                   *
   6  *                                                                        *
   7  * Primary Authors: Oystein Djuvsland                                     *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 /**
  19  * @file   AliHLTPHOSClusterizer.cxx
  20  * @author Oystein Djuvsland
  21  * @date
  22  * @brief  Clusterizer for PHOS HLT
  23  */
  24
  25 // see header file for class documentation
  26 // or
  27 // refer to README to build package
  28 // or
  29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  30
  31 #include "AliHLTPHOSClusterizer.h"
  32 //#include "AliHLTPHOSBase.h"
  33 #include "AliHLTLogging.h"
  34 #include "TMath.h"
  35 #include "AliHLTPHOSRecPointContainerStruct.h"
  36 #include "AliHLTPHOSRecPointDataStruct.h"
  37 #include "AliHLTPHOSDigitDataStruct.h"
  38 #include "AliHLTPHOSDigitContainerDataStruct.h"
  39 #include "AliHLTPHOSDigitReader.h"
  40 #include "TClonesArray.h"
  41 #include "AliPHOSDigit.h"
  42 #ifndef HAVENOT__PHOSRECOPARAMEMC // set from configure if EMC functionality not available in AliPHOSRecoParam
  43 #include "AliPHOSRecoParam.h"
  44 #else
  45 #include "AliPHOSRecoParamEmc.h"
  46 #endif
  47 #include <iostream>
  48
  49 using namespace std;
  50
  51 ClassImp(AliHLTPHOSClusterizer);
  52
  53 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
  54   AliHLTLogging(),
  55   fRecPointDataPtr(0),
  56   fDigitDataPtr(0),
  57   fCurrentDigit(0),
  58   fStartDigit(0),
  59   fPreviousDigit(0),
  60   fEmcClusteringThreshold(0),
  61   fEmcMinEnergyThreshold(0),
  62   fEmcTimeGate(0),
  63   fDigitsInCluster(0),
  64   fDigitContainerPtr(0),
  65   fMaxDigitIndexDiff(2*NZROWSMOD),
  66   fAvailableSize(0),
  67   fDigitReader(0)
  68 {
  69   //See header file for documentation
  70   fEmcClusteringThreshold = 0.2;
  71   fEmcMinEnergyThreshold = 0.03;
  72   fEmcTimeGate = 1.e-6 ;
  73
  74   fDigitReader = new AliHLTPHOSDigitReader();
  75 }//end
  76
  77
  78 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()
  79 {
  80   //See header file for documentation
  81 }
  82
  83 void
  84 AliHLTPHOSClusterizer::SetRecPointDataPtr(AliHLTPHOSRecPointDataStruct* recPointDataPtr)
  85 {
  86   // See header file for documentation
  87   fRecPointDataPtr = recPointDataPtr;
  88 }
  89
  90 Int_t
  91 AliHLTPHOSClusterizer::ClusterizeEvent(AliHLTPHOSDigitHeaderStruct *digitHeader, UInt_t availableSize, UInt_t& totSize)
  92 {
  93   //see header file for documentation
  94   Int_t nRecPoints = 0;
  95
  96   fAvailableSize = availableSize;
  97
  98   fDigitReader->SetDigitHeader(digitHeader);
  99
 100   //  UInt_t maxRecPointSize = sizeof(AliHLTPHOSRecPointDataStruct) + (sizeof(AliHLTPHOSDigitDataStruct) << 7); //Reasonable estimate...
 101
 102   //  HLTError("Starting clusterisation");
 103
 104   //Clusterization starts
 105   while((fCurrentDigit = fDigitReader->NextDigit()) != 0)
 106     {
 107       //      HLTError("Digit with energy: %f", fCurrentDigit->fEnergy);
 108       fDigitsInCluster = 0;
 109
 110       if(fCurrentDigit->fEnergy < fEmcClusteringThreshold)
 111         {
 112           continue;
 113         }
 114
 115       //      HLTError("Have cluster candidate (x,z): %d, %d, with energy: %f", fCurrentDigit->fX, fCurrentDigit->fZ, fCurrentDigit->fEnergy);
 116
 117       if(fAvailableSize < (sizeof(AliHLTPHOSRecPointDataStruct)))
 118         {
 119           HLTError("Out of buffer, stopping clusterisation");
 120           return -1;
 121         }
 122
 123       // Save the starting digit
 124       fStartDigit = fCurrentDigit;
 125       fPreviousDigit = fStartDigit;
 126       // Get the offset of the digit starting the cluster relative to the start of the digit block
 127       fRecPointDataPtr->fStartDigitOffset = fDigitReader->GetCurrentDigitOffset();
 128       //      HLTError("Start digit offset: %d", fRecPointDataPtr->fStartDigitOffset);
 129       //      cout << "Start digit offset: " <<  fRecPointDataPtr->fStartDigitOffset << endl;
 130       fDigitReader->DropDigit();
 131
 132       fRecPointDataPtr->fAmp = 0;
 133       fRecPointDataPtr->fModule = fCurrentDigit->fModule;
 134
 135       fRecPointDataPtr->fAmp += fCurrentDigit->fEnergy;
 136       fDigitsInCluster++;
 137
 138       nRecPoints++;
 139
 140       // Scanning for the neighbours
 141       if(ScanForNeighbourDigits(fRecPointDataPtr, fCurrentDigit) < 0)
 142         {
 143           return -1;
 144         }
 145
 146       fPreviousDigit->fMemOffsetNext = 0;
 147
 148       totSize += sizeof(AliHLTPHOSRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTPHOSDigitDataStruct);
 149       HLTDebug("Initial available size: %d, used size: %d, remaining available size: %d, should be: %d", availableSize, totSize, fAvailableSize, availableSize-totSize);
 150
 151       fRecPointDataPtr->fMultiplicity = fDigitsInCluster;
 152       HLTDebug("Number of digits in cluster: %d", fDigitsInCluster);
 153       fRecPointDataPtr++;
 154
 155       fDigitReader->Rewind(); // TODO: jump to the next instead of rewind
 156     }//end of clusterization
 157
 158   return nRecPoints;
 159 }
 160
 161 Int_t
 162 AliHLTPHOSClusterizer::ScanForNeighbourDigits(AliHLTPHOSRecPointDataStruct* recPoint, AliHLTPHOSDigitDataStruct *digit)
 163 {
 164   //see header file for documentation
 165
 166   //Int_t max = TMath::Min((Int_t)fDigitContainerPtr->fNDigits, (Int_t)fMaxDigitIndexDiff+index);
 167   //  Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));
 168
 169   //  Int_t max = fDigitContainerPtr->fNDigits;
 170   //  Int_t min = 0;
 171
 172   AliHLTPHOSDigitDataStruct *tmpDigit = 0;
 173
 174   fDigitReader->Rewind();
 175
 176   while((tmpDigit = fDigitReader->NextDigit()))
 177     {
 178       //HLTError("Checking digit (x,z): %d, %d, with energy: %f", tmpDigit->fX, tmpDigit->fZ, tmpDigit->fEnergy);
 179       //      HLTError("Checking digit (x,z): %d, %d, with energy: %f for neighbourship of digit (x,z): %d, %d", tmpDigit->fX, tmpDigit->fZ, tmpDigit->fEnergy, digit->fX, digit->fZ);
 180       //      if(tmpDigit->fEnergy > fEmcMinEnergyThreshold)
 181       if(tmpDigit->fEnergy > 0.2)
 182         {
 183           //HLTError("Checking digit (x,z): %d, %d, with energy: %f for neighbourship of digit (x,z): %d, %d", tmpDigit->fX, tmpDigit->fZ, tmpDigit->fEnergy, digit->fX, digit->fZ);
 184
 185           //      if(AreNeighbours(fStartDigit, tmpDigit))
 186           if(AreNeighbours(digit, tmpDigit))
 187             {
 188               //              HLTError("Adding digit (x,z): %d, %d, with energy: %f", tmpDigit->fX, tmpDigit->fZ, tmpDigit->fEnergy);
 189               fDigitReader->DropDigit();
 190               fPreviousDigit->fMemOffsetNext = reinterpret_cast<Long_t>(tmpDigit) - reinterpret_cast<Long_t>(fPreviousDigit);
 191               //              cout << "Digit offset: " << fPreviousDigit->fMemOffsetNext << endl;
 192               //   fPreviousDigit->fMemOffsetNext += tmpDigit->fMemOffsetNext;
 193               recPoint->fAmp += tmpDigit->fEnergy;
 194               fPreviousDigit = tmpDigit;
 195               fDigitsInCluster++;
 196               ScanForNeighbourDigits(recPoint, tmpDigit);
 197               fDigitReader->SetCurrentDigit(tmpDigit);
 198             }
 199         }
 200     }
 201   return 0;
 202 }
 203
 204 Int_t
 205 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1,
 206                                             AliHLTPHOSDigitDataStruct* digit2)
 207 {
 208   //see header file for documentation
 209   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
 210     {
 211       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
 212       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
 213       //      HLTError("coldiff: %d, rowdiff: %d, timediff: %f", coldiff, rowdiff, TMath::Abs(digit1->fTime - digit2->fTime ));
 214       if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))
 215         {
 216           //      cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ <<
 217             //      " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;
 218
 219           //      if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
 220             {
 221               return 1;
 222             }
 223         }
 224     }
 225   return 0;
 226 }