HLT/PHOS/AliHLTPHOSClusterizer.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * All rights reserved.                                                   *
   4  *                                                                        *
   5  * Primary Authors: Oystein Djuvsland                                     *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /**
  17  * @file   AliHLTPHOSClusterizer.cxx
  18  * @author Oystein Djuvsland
  19  * @date
  20  * @brief  Clusterizer for PHOS HLT
  21  */
  22
  23 // see header file for class documentation
  24 // or
  25 // refer to README to build package
  26 // or
  27 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  28
  29 #include "AliHLTPHOSClusterizer.h"
  30 #include "AliHLTPHOSBase.h"
  31 #include "TMath.h"
  32 #include "AliHLTPHOSRecPointContainerStruct.h"
  33 #include "AliHLTPHOSRecPointDataStruct.h"
  34 #include "AliHLTPHOSDigitDataStruct.h"
  35 #include "AliHLTPHOSDigitContainerDataStruct.h"
  36 #include "TClonesArray.h"
  37 #include "AliPHOSGeometry.h"
  38 #include "AliPHOSDigit.h"
  39 #include "AliPHOSRecoParam.h"
  40
  41 ClassImp(AliHLTPHOSClusterizer);
  42
  43 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
  44   AliHLTPHOSBase(),
  45   fEmcClusteringThreshold(0),
  46   fEmcMinEnergyThreshold(0),
  47   fEmcTimeGate(0),
  48   fLogWeight(0),
  49   fDigitsInCluster(0),
  50   fOnlineMode(true),
  51   fDigitArrayPtr(0),
  52   fEmcRecPointsPtr(0),
  53   fDigitPtr(0),
  54   fDigitContainerPtr(0),
  55   fRecPointContainerPtr(0),
  56   fPHOSGeometry(0),
  57   fGetterPtr(0)
  58 {
  59   //See header file for documentation
  60   fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
  61   fEmcClusteringThreshold = 0.2;
  62   fEmcMinEnergyThreshold = 0.03;
  63   fEmcTimeGate = 1.e-6 ;
  64   fLogWeight = 4.5;
  65
  66
  67 }//end
  68
  69
  70 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()
  71 {
  72   //See header file for documentation
  73 }
  74
  75 void
  76 AliHLTPHOSClusterizer::SetRecPointContainer(AliHLTPHOSRecPointContainerStruct* recPointContainerPtr)
  77   {
  78     fRecPointContainerPtr = recPointContainerPtr;
  79     fRecPointContainerPtr->fNRecPoints = 0;
  80   }
  81
  82 void
  83 AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParam* params)
  84 {
  85   //see header file for documentation
  86   fEmcClusteringThreshold = params->GetEMCClusteringThreshold();
  87   fEmcMinEnergyThreshold = params->GetEMCMinE();
  88   fLogWeight = params->GetEMCLogWeight();
  89 }
  90
  91 void
  92 AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSLoader* getter)
  93 {
  94   //see header file for documentation
  95   fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
  96   fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
  97   fGetterPtr = getter;
  98   fDigitArrayPtr = fGetterPtr->Digits();
  99   fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
 100   fOnlineMode = false;
 101 }
 102
 103 Int_t
 104 AliHLTPHOSClusterizer::GetEvent(Int_t i)
 105 {
 106   //see header file for documentation
 107   Int_t coord[4];
 108
 109   //  fGetterPtr->Event(i, "D");
 110   fGetterPtr->GetEvent();
 111   for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
 112     {
 113       fDigitPtr = (AliPHOSDigit*)fDigitArrayPtr->At(j);
 114       fPHOSGeometry->AbsToRelNumbering(fDigitPtr->GetId(), coord);
 115       fDigitContainerPtr->fDigitDataStruct[j].fX = coord[3];
 116       fDigitContainerPtr->fDigitDataStruct[j].fZ = coord[2];
 117       fDigitContainerPtr->fDigitDataStruct[j].fModule = coord[0];
 118       fDigitContainerPtr->fDigitDataStruct[j].fEnergy = fDigitPtr->GetEnergy();
 119       fDigitContainerPtr->fDigitDataStruct[j].fTime = fDigitPtr->GetTime();
 120     }
 121   fDigitContainerPtr->fNDigits = fDigitArrayPtr->GetEntriesFast();
 122   return 0;
 123 }
 124
 125 Int_t
 126 AliHLTPHOSClusterizer::GetNEvents()
 127 {
 128   //see header file for documentation
 129   if(fOnlineMode)
 130     {
 131       //     Logging(kHLTLogWarning, __FILE__ , "information not available" , "GetNEvents()  Number of events not available in online mod");
 132       return -1;
 133     }
 134   return fGetterPtr->MaxEvent();
 135 }
 136
 137
 138 Int_t
 139 AliHLTPHOSClusterizer::ClusterizeEvent()
 140 {
 141   //see header file for documentation
 142   Int_t nRecPoints = 0;
 143   UInt_t i = 0;
 144
 145   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 146
 147   //Clusterization starts
 148   for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
 149     {
 150
 151       fDigitsInCluster = 0;
 152       if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
 153         {
 154           continue;
 155         }
 156       recPoint = &(fRecPointContainerPtr->fRecPointArray[nRecPoints]);
 157       recPoint->fAmp = 0;
 158       //TODO!!!!!!!
 159       recPoint->fModule = fDigitContainerPtr->fDigitDataStruct[i].fModule;
 160       //TODO!!!!!!!
 161       //recPoint->
 162       recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[i];
 163       recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[i].fEnergy;
 164       fDigitContainerPtr->fDigitDataStruct[i].fEnergy = 0;
 165       fDigitsInCluster++;
 166       nRecPoints++;
 167       ScanForNeighbourDigits(i, recPoint);
 168
 169       recPoint->fMultiplicity = fDigitsInCluster;
 170
 171     }//end of clusterization
 172   fRecPointContainerPtr->fNRecPoints = nRecPoints;
 173
 174   return nRecPoints;
 175 }
 176
 177 void
 178 AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
 179
 180 {
 181   //see header file for documentation
 182
 183   for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
 184     {
 185       if(fDigitContainerPtr->fDigitDataStruct[j].fEnergy > fEmcMinEnergyThreshold)
 186         {
 187           switch(AreNeighbours(&(fDigitContainerPtr->fDigitDataStruct[index]),
 188                                &(fDigitContainerPtr->fDigitDataStruct[j])))
 189             {
 190             case 0:
 191               break;
 192             case 1:
 193               recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[j];
 194               recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[j].fEnergy;
 195               fDigitContainerPtr->fDigitDataStruct[j].fEnergy = 0;
 196               fDigitsInCluster++;
 197               ScanForNeighbourDigits(j, recPoint);
 198               break;
 199             case 2:
 200               break;
 201             }
 202         }
 203     }
 204   return;
 205 }
 206
 207
 208 Int_t
 209 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1,
 210                                             AliHLTPHOSDigitDataStruct* digit2)
 211 {
 212   //see header file for documentation
 213
 214   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
 215     {
 216       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
 217       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
 218
 219       if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
 220         {
 221           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
 222             {
 223               return 1;
 224             }
 225         }
 226     }
 227   return 0;
 228 }
 229
 230 void
 231 AliHLTPHOSClusterizer::CalculateCenterOfGravity()
 232 {
 233   //see header file for documentation
 234   Float_t wtot = 0.;
 235   Float_t x = 0.;
 236   Float_t z = 0.;
 237   Float_t xi = 0.;
 238   Float_t zi = 0.;
 239
 240   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 241   AliHLTPHOSDigitDataStruct *digit = 0;
 242   UInt_t iDigit = 0;
 243   UInt_t iRecPoint = 0;
 244
 245   for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++)
 246     {
 247       recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
 248       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 249         {
 250           digit = &(recPoint->fDigitsList[iDigit]);
 251
 252           //fPHOSGeometry->AbsToRelNumbering(digit->fID, relid) ;
 253           //  fPHOSGeometry->RelPosInModule(relid, xi, zi);
 254           xi = digit->fX;
 255           zi = digit->fZ;
 256
 257           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 258             {
 259               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 260               x    += xi * w ;
 261               z    += zi * w ;
 262               wtot += w ;
 263             }
 264         }
 265
 266       if (wtot>0)
 267         {
 268           recPoint->fX = x/wtot ;
 269           recPoint->fZ = z/wtot ;
 270         }
 271       else
 272         {
 273           recPoint->fAmp = 0;
 274         }
 275     }
 276
 277 }
 278
 279