HLT/PHOS/AliHLTPHOSClusterizer.cxx

   1 //Insert copyright
   2
   3
   4 #include "AliHLTPHOSClusterizer.h"
   5 #include "AliHLTPHOSBase.h"
   6 #include "TMath.h"
   7 #include "AliHLTPHOSRecPointContainerStruct.h"
   8 #include "AliHLTPHOSRecPointDataStruct.h"
   9 #include "AliHLTPHOSDigitDataStruct.h"
  10 #include "AliHLTPHOSDigitContainerDataStruct.h"
  11 #include "TClonesArray.h"
  12 #include "AliPHOSGeometry.h"
  13 #include "AliPHOSDigit.h"
  14 #include "AliPHOSRecoParamEmc.h"
  15
  16 ClassImp(AliHLTPHOSClusterizer);
  17
  18 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
  19   AliHLTPHOSBase(),
  20   fEmcClusteringThreshold(0),
  21   fEmcMinEnergyThreshold(0),
  22   fEmcTimeGate(0),
  23   fLogWeight(0),
  24   fDigitsInCluster(0),
  25   fOnlineMode(true),
  26   fDigitArrayPtr(0),
  27   fEmcRecPointsPtr(0),
  28   fDigitPtr(0),
  29   fDigitContainerPtr(0),
  30   fRecPointContainerPtr(0),
  31   fPHOSGeometry(0),
  32   fGetterPtr(0)
  33 {
  34   //See header file for documentation
  35   fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
  36   fEmcClusteringThreshold = 0.2;
  37   fEmcMinEnergyThreshold = 0.03;
  38   fEmcTimeGate = 1.e-8 ;
  39   fLogWeight = 4.5;
  40
  41
  42 }//end
  43
  44
  45 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &)
  46 {
  47   //Copy constructor, not implemented
  48 }//end
  49
  50 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()
  51 {
  52   //Destructor
  53 }
  54
  55 void
  56 AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParamEmc* params)
  57 {
  58   fEmcClusteringThreshold = params->GetClusteringThreshold();
  59   fEmcMinEnergyThreshold = params->GetMinE();
  60   fLogWeight = params->GetLogWeight();
  61 }
  62
  63 void
  64 AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSGetter* getter)
  65 {
  66   fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
  67   fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
  68   fGetterPtr = getter;
  69   fDigitArrayPtr = fGetterPtr->Digits();
  70   fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
  71   fOnlineMode = false;
  72 }
  73
  74 Int_t
  75 AliHLTPHOSClusterizer::GetEvent(Int_t i)
  76 {
  77   Int_t coord[4];
  78
  79   fGetterPtr->Event(i, "D");
  80   for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
  81     {
  82       fDigitPtr = (AliPHOSDigit*)fDigitArrayPtr->At(j);
  83       fPHOSGeometry->AbsToRelNumbering(fDigitPtr->GetId(), coord);
  84       fDigitContainerPtr->fDigitDataStruct[j].fX = coord[3];
  85       fDigitContainerPtr->fDigitDataStruct[j].fZ = coord[2];
  86       fDigitContainerPtr->fDigitDataStruct[j].fModule = coord[0];
  87       fDigitContainerPtr->fDigitDataStruct[j].fEnergy = fDigitPtr->GetEnergy();
  88       fDigitContainerPtr->fDigitDataStruct[j].fTime = fDigitPtr->GetTime();
  89     }
  90   fDigitContainerPtr->fNDigits = fDigitArrayPtr->GetEntriesFast();
  91   return 0;
  92 }
  93
  94 Int_t
  95 AliHLTPHOSClusterizer::GetNEvents()
  96 {
  97   if(fOnlineMode)
  98     {
  99       printf("Number of events not available in online mode!\n");
 100       return -1;
 101     }
 102   return fGetterPtr->MaxEvent();
 103 }
 104
 105
 106 Int_t
 107 AliHLTPHOSClusterizer::ClusterizeEvent()
 108 {
 109   // Steering method to construct the clusters stored in a list of Reconstructed Points
 110   // A cluster is defined as a list of neighbour digits
 111   Int_t nRecPoints = 0;
 112   UInt_t i = 0;
 113
 114
 115   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 116
 117   //Clusterization starts
 118   for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
 119     {
 120
 121       fDigitsInCluster = 0;
 122
 123       if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
 124         {
 125           continue;
 126         }
 127
 128       recPoint = &(fRecPointContainerPtr->fRecPointArray[nRecPoints]);
 129       recPoint->fAmp = 0;
 130       //recPoint->
 131       recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[i];
 132       recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[i].fEnergy;
 133       fDigitContainerPtr->fDigitDataStruct[i].fEnergy = 0;
 134       fDigitsInCluster++;
 135       nRecPoints++;
 136       ScanForNeighbourDigits(i, recPoint);
 137
 138       recPoint->fMultiplicity = fDigitsInCluster;
 139     }//end of clusterization
 140   fRecPointContainerPtr->fNRecPoints = nRecPoints;
 141
 142   return nRecPoints;
 143 }
 144
 145 void
 146 AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
 147
 148 {
 149
 150   for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
 151     {
 152       if(fDigitContainerPtr->fDigitDataStruct[j].fEnergy > fEmcMinEnergyThreshold)
 153         {
 154           switch(AreNeighbours(&(fDigitContainerPtr->fDigitDataStruct[index]),
 155                                &(fDigitContainerPtr->fDigitDataStruct[j])))
 156             {
 157             case 0:
 158               break;
 159             case 1:
 160               recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[j];
 161               recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[j].fEnergy;
 162               fDigitContainerPtr->fDigitDataStruct[j].fEnergy = 0;
 163               fDigitsInCluster++;
 164               ScanForNeighbourDigits(j, recPoint);
 165               break;
 166             case 2:
 167               break;
 168             }
 169         }
 170     }
 171   return;
 172 }
 173
 174
 175 Int_t
 176 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1,
 177                                             AliHLTPHOSDigitDataStruct* digit2)
 178 {
 179
 180   //Int_t coord1[4];
 181   //Int_t coord2[4];
 182
 183   // fPHOSGeometry->AbsToRelNumbering(digit1->fID, coord1);
 184   //fPHOSGeometry->AbsToRelNumbering(digit2->fID, coord2);
 185
 186
 187   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
 188     {
 189
 190       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
 191       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
 192
 193       if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
 194         {
 195           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
 196             return 1;
 197         }
 198     }
 199   return 0;
 200 }
 201
 202 void
 203 AliHLTPHOSClusterizer::CalculateCenterOfGravity()
 204 {
 205
 206   // Calculates the center of gravity in the local PHOS-module coordinates
 207   Float_t wtot = 0.;
 208
 209   Int_t relid[4];
 210
 211   Float_t x = 0.;
 212   Float_t z = 0.;
 213   Float_t xi = 0.;
 214   Float_t zi = 0.;
 215
 216   AliHLTPHOSRecPointDataStruct *recPoint = 0;
 217   AliHLTPHOSDigitDataStruct *digit = 0;
 218
 219   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 220
 221   Int_t iDigit = 0;
 222   Int_t iRecPoint = 0;
 223
 224   for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++)
 225     {
 226       recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
 227       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 228         {
 229           digit = &(recPoint->fDigitsList[iDigit]);
 230
 231           //fPHOSGeometry->AbsToRelNumbering(digit->fID, relid) ;
 232           //  fPHOSGeometry->RelPosInModule(relid, xi, zi);
 233           xi = digit->fX;
 234           zi = digit->fZ;
 235
 236           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 237             {
 238               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 239               x    += xi * w ;
 240               z    += zi * w ;
 241               wtot += w ;
 242             }
 243         }
 244
 245       if (wtot>0)
 246         {
 247           recPoint->fX = x/wtot ;
 248           recPoint->fZ = z/wtot ;
 249         }
 250       else
 251         {
 252           recPoint->fAmp = 0;
 253         }
 254     }
 255
 256 }
 257
 258