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 #ifndef HAVE_NOT_PHOSRECOPARAMEMC // set from configure if EMC functionality not available in AliPHOSRecoParam
  40 #include "AliPHOSRecoParam.h"
  41 #else
  42 #include "AliPHOSRecoParamEmc.h"
  43 #endif
  44
  45 ClassImp(AliHLTPHOSClusterizer);
  46
  47 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
  48   AliHLTPHOSBase(),
  49   fEmcClusteringThreshold(0),
  50   fEmcMinEnergyThreshold(0),
  51   fEmcTimeGate(0),
  52   fLogWeight(0),
  53   fDigitsInCluster(0),
  54   fOnlineMode(true),
  55   fDigitArrayPtr(0),
  56   fEmcRecPointsPtr(0),
  57   fDigitPtr(0),
  58   fDigitContainerPtr(0),
  59   fRecPointContainerPtr(0),
  60   fPHOSGeometry(0),
  61   fGetterPtr(0)
  62 {
  63   //See header file for documentation
  64   fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
  65   fEmcClusteringThreshold = 0.2;
  66   fEmcMinEnergyThreshold = 0.03;
  67   fEmcTimeGate = 1.e-6 ;
  68   fLogWeight = 4.5;
  69
  70
  71 }//end
  72
  73
  74 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()
  75 {
  76   //See header file for documentation
  77 }
  78
  79 void
  80 AliHLTPHOSClusterizer::SetRecPointContainer(AliHLTPHOSRecPointContainerStruct* recPointContainerPtr)
  81   {
  82     fRecPointContainerPtr = recPointContainerPtr;
  83     fRecPointContainerPtr->fNRecPoints = 0;
  84   }
  85
  86 void
  87 AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParam* params)
  88 {
  89   //see header file for documentation
  90 #ifndef HAVE_NOT_PHOSRECOPARAMEMC // set from configure if EMC functionality not available in AliPHOSRecoParam
  91   // the new AliPHOSRecoParam functions, available from revision
  92   fEmcClusteringThreshold = params->GetEMCClusteringThreshold();
  93   fEmcMinEnergyThreshold = params->GetEMCMinE();
  94   fLogWeight = params->GetEMCLogWeight();
  95 #else
  96   fEmcClusteringThreshold = params->GetClusteringThreshold();
  97   fEmcMinEnergyThreshold = params->GetMinE();
  98   fLogWeight = params->GetLogWeight();
  99 #endif
 100 }
 101
 102 void
 103 AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSLoader* getter)
 104 {
 105   //see header file for documentation
 106   fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
 107   fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
 108   fGetterPtr = getter;
 109   fDigitArrayPtr = fGetterPtr->Digits();
 110   fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
 111   fOnlineMode = false;
 112 }
 113
 114 Int_t
 115 AliHLTPHOSClusterizer::GetEvent(Int_t i)
 116 {
 117   //see header file for documentation
 118   Int_t coord[4];
 119
 120   //  fGetterPtr->Event(i, "D");
 121   fGetterPtr->GetEvent();
 122   for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
 123     {
 124       fDigitPtr = (AliPHOSDigit*)fDigitArrayPtr->At(j);
 125       fPHOSGeometry->AbsToRelNumbering(fDigitPtr->GetId(), coord);
 126       fDigitContainerPtr->fDigitDataStruct[j].fX = coord[3];
 127       fDigitContainerPtr->fDigitDataStruct[j].fZ = coord[2];
 128       fDigitContainerPtr->fDigitDataStruct[j].fModule = coord[0];
 129       fDigitContainerPtr->fDigitDataStruct[j].fEnergy = fDigitPtr->GetEnergy();
 130       fDigitContainerPtr->fDigitDataStruct[j].fTime = fDigitPtr->GetTime();
 131     }
 132   fDigitContainerPtr->fNDigits = fDigitArrayPtr->GetEntriesFast();
 133   return 0;
 134 }
 135
 136 Int_t
 137 AliHLTPHOSClusterizer::GetNEvents()
 138 {
 139   //see header file for documentation
 140   if(fOnlineMode)
 141     {
 142       //     Logging(kHLTLogWarning, __FILE__ , "information not available" , "GetNEvents()  Number of events not available in online mod");
 143       return -1;
 144     }
 145   return fGetterPtr->MaxEvent();
 146 }
 147
 148
 149 Int_t
 150 AliHLTPHOSClusterizer::ClusterizeEvent()
 151 {
 152   //see header file for documentation
 153   Int_t nRecPoints = 0;
 154   UInt_t i = 0;
 155
 156   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 157
 158   //Clusterization starts
 159   for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
 160     {
 161
 162       fDigitsInCluster = 0;
 163       if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
 164         {
 165           continue;
 166         }
 167       recPoint = &(fRecPointContainerPtr->fRecPointArray[nRecPoints]);
 168       recPoint->fAmp = 0;
 169       //TODO!!!!!!!
 170       recPoint->fModule = fDigitContainerPtr->fDigitDataStruct[i].fModule;
 171       //TODO!!!!!!!
 172       //recPoint->
 173       recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[i];
 174       recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[i].fEnergy;
 175       fDigitContainerPtr->fDigitDataStruct[i].fEnergy = 0;
 176       fDigitsInCluster++;
 177       nRecPoints++;
 178       ScanForNeighbourDigits(i, recPoint);
 179
 180       recPoint->fMultiplicity = fDigitsInCluster;
 181
 182     }//end of clusterization
 183   fRecPointContainerPtr->fNRecPoints = nRecPoints;
 184
 185   return nRecPoints;
 186 }
 187
 188 void
 189 AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
 190
 191 {
 192   //see header file for documentation
 193
 194   for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
 195     {
 196       if(fDigitContainerPtr->fDigitDataStruct[j].fEnergy > fEmcMinEnergyThreshold)
 197         {
 198           switch(AreNeighbours(&(fDigitContainerPtr->fDigitDataStruct[index]),
 199                                &(fDigitContainerPtr->fDigitDataStruct[j])))
 200             {
 201             case 0:
 202               break;
 203             case 1:
 204               recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[j];
 205               recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[j].fEnergy;
 206               fDigitContainerPtr->fDigitDataStruct[j].fEnergy = 0;
 207               fDigitsInCluster++;
 208               ScanForNeighbourDigits(j, recPoint);
 209               break;
 210             case 2:
 211               break;
 212             }
 213         }
 214     }
 215   return;
 216 }
 217
 218
 219 Int_t
 220 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1,
 221                                             AliHLTPHOSDigitDataStruct* digit2)
 222 {
 223   //see header file for documentation
 224
 225   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
 226     {
 227       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
 228       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
 229
 230       if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
 231         {
 232           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
 233             {
 234               return 1;
 235             }
 236         }
 237     }
 238   return 0;
 239 }
 240
 241 void
 242 AliHLTPHOSClusterizer::CalculateCenterOfGravity()
 243 {
 244   //see header file for documentation
 245   Float_t wtot = 0.;
 246   Float_t x = 0.;
 247   Float_t z = 0.;
 248   Float_t xi = 0.;
 249   Float_t zi = 0.;
 250
 251   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 252   AliHLTPHOSDigitDataStruct *digit = 0;
 253   UInt_t iDigit = 0;
 254   UInt_t iRecPoint = 0;
 255
 256   for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++)
 257     {
 258       recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
 259       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 260         {
 261           digit = &(recPoint->fDigitsList[iDigit]);
 262
 263           //fPHOSGeometry->AbsToRelNumbering(digit->fID, relid) ;
 264           //  fPHOSGeometry->RelPosInModule(relid, xi, zi);
 265           xi = digit->fX;
 266           zi = digit->fZ;
 267
 268           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 269             {
 270               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 271               x    += xi * w ;
 272               z    += zi * w ;
 273               wtot += w ;
 274             }
 275         }
 276
 277       if (wtot>0)
 278         {
 279           recPoint->fX = x/wtot ;
 280           recPoint->fZ = z/wtot ;
 281         }
 282       else
 283         {
 284           recPoint->fAmp = 0;
 285         }
 286     }
 287
 288 }
 289
 290