HLT/CALO/AliHLTCaloClusterAnalyser.cxx

   1 // $Id: AliHLTCaloClusterAnalyser.cxx 35107 2009-09-30 01:45:06Z phille $
   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   AliHLTCaloClusterAnalyser.cxx
  20  * @author Oystein Djuvsland
  21  * @date
  22  * @brief  Cluster analyser for Calo 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 "AliHLTCaloClusterAnalyser.h"
  32 #include "AliHLTCaloRecPointHeaderStruct.h"
  33 #include "AliHLTCaloRecPointDataStruct.h"
  34 #include "AliHLTCaloClusterDataStruct.h"
  35 //#include "AliHLTCaloPhysicsAnalyzer.h"
  36 #include "AliHLTCaloGeometry.h"
  37 #include "AliESDCaloCluster.h"
  38 #include "TMath.h"
  39 #include "TVector3.h"
  40 #include "TH1F.h"
  41 #include "TFile.h"
  42 #include "AliHLTCaloRecoParamHandler.h"
  43
  44 ClassImp(AliHLTCaloClusterAnalyser);
  45
  46 AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
  47   //  AliHLTCaloBase(),
  48   fLogWeight(4.5),
  49   fRecPointArray(0),
  50   fDigitDataArray(0),
  51   fNRecPoints(0),
  52   fCaloClusterDataPtr(0),
  53   fCaloClusterHeaderPtr(0),
  54   //fAnalyzerPtr(0),
  55   fDoClusterFit(false),
  56   fHaveCPVInfo(false),
  57   fDoPID(false),
  58   fHaveDistanceToBadChannel(false),
  59   fGeometry(0),
  60   fClusterType(AliESDCaloCluster::kPHOSCluster),
  61   fRecoParamsPtr(0)
  62 {
  63   //See header file for documentation
  64 }
  65
  66 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
  67 {
  68   // See header file for class documentation
  69 }
  70
  71 void
  72 AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
  73 {
  74   //see header file for documentation
  75   fCaloClusterDataPtr = caloClusterDataPtr;
  76 }
  77
  78 void
  79 AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
  80 {
  81   fRecPointArray = recPointDataPtr;
  82   fNRecPoints = nRecPoints;
  83 }
  84
  85 void
  86 AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits)
  87 {
  88 //   AliHLTCaloClusterizer cl("PHOS");
  89   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
  90    fDigitDataArray = digits;
  91    //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
  92 }
  93
  94 Int_t
  95 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
  96 {
  97   //see header file for documentation
  98   Float_t wtot = 0.;
  99   Float_t x = 0.;
 100   Float_t z = 0.;
 101   Float_t xi = 0.;
 102   Float_t zi = 0.;
 103
 104   AliHLTCaloDigitDataStruct *digit = 0;
 105
 106   UInt_t iDigit = 0;
 107
 108   for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
 109     {
 110       AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
 111       //      digit = &(recPoint->fDigits);
 112
 113       Int_t *digitIndexPtr = &(recPoint->fDigits);
 114        wtot = 0;
 115        x = 0;
 116        z = 0;
 117
 118 /*       Float_t maxAmp = 0;
 119        Int_t maxX = 0;
 120        Int_t maxZ = 0;*/
 121       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 122         {
 123
 124           digit = &(fDigitDataArray[*digitIndexPtr]);
 125
 126           xi = digit->fX+0.5;
 127           zi = digit->fZ+0.5;
 128
 129           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 130             {
 131               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 132               x    += xi * w ;
 133               z    += zi * w ;
 134               wtot += w ;
 135 /*            if(digit->fEnergy > maxAmp)
 136               {
 137                  maxAmp = digit->fEnergy;
 138                  maxX = digit->fX;// + 0.5;
 139                  maxZ = digit->fZ;// + 0.5;
 140               }*/
 141             }
 142           digitIndexPtr++;
 143         }
 144
 145       if (wtot>0)
 146         {
 147           recPoint->fX = x/wtot ;
 148           recPoint->fZ = z/wtot ;
 149         }
 150       else
 151         {
 152           recPoint->fAmp = 0;
 153         }
 154 //      printf("Max digit: E = %f, x = %d, z= %d, cluster: E = %f, x = %f, z = %f\n" , maxAmp, maxX, maxZ, recPoint->fAmp, recPoint->fX, recPoint->fZ);
 155     }
 156   return 0;
 157 }
 158
 159
 160 Int_t
 161 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
 162 {
 163   //See header file for documentation
 164   return 0;
 165 }
 166
 167 Int_t
 168 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
 169 {
 170   //See header file for documentation
 171   return 0;
 172 }
 173
 174
 175 Int_t
 176 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
 177 {
 178   //See header file for documentation
 179   return 0;
 180 }
 181
 182 Int_t
 183 AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
 184 {
 185   //See header file for documentation
 186
 187    if(fRecoParamsPtr)
 188    {
 189       fLogWeight = fRecoParamsPtr->GetLogWeight();
 190    }
 191
 192   fNRecPoints = nRecPoints;
 193
 194   if(fGeometry == 0)
 195   {
 196      HLTError("No geometry object is initialised, creation of clusters stopped");
 197   }
 198
 199   CalculateCenterOfGravity();
 200
 201   //  AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
 202   AliHLTCaloDigitDataStruct* digitPtr = 0;
 203
 204   AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
 205
 206   //  Int_t id = -1;
 207   TVector3 globalPos;
 208
 209   for(Int_t i = 0; i < fNRecPoints; i++) //TODO needs fix when we start unfolding (number of clusters not necessarily same as number of recpoints gotten from the clusterizer
 210     {
 211       if((availableSize - totSize)  < sizeof(AliHLTCaloClusterDataStruct))
 212       {
 213          HLTError("Out of buffer: available size is: %d, total size used: %d", availableSize, totSize);
 214          return -ENOBUFS;
 215       }
 216       totSize += sizeof(AliHLTCaloClusterDataStruct);
 217
 218       caloClusterPtr = fCaloClusterDataPtr;
 219       caloClusterPtr->fChi2 = 0;
 220       caloClusterPtr->fClusterType = kUndef;
 221       caloClusterPtr->fDispersion = 0;
 222       caloClusterPtr->fDistanceToBadChannel = 0;
 223       caloClusterPtr->fDistToBadChannel = 0;
 224       caloClusterPtr->fEmcCpvDistance = 0;
 225       caloClusterPtr->fEnergy = 0;
 226       caloClusterPtr->fFitQuality = 0;
 227       caloClusterPtr->fID = 0;
 228       caloClusterPtr->fM02 = 0;
 229       caloClusterPtr->fM20 = 0;
 230       caloClusterPtr->fNCells = 0;
 231       caloClusterPtr->fNExMax = 0;
 232       caloClusterPtr->fTOF = 0;
 233       caloClusterPtr->fTrackDx = 0;
 234       caloClusterPtr->fTrackDz = 0;
 235
 236       AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
 237
 238       AliHLTCaloGlobalCoordinate globalCoord;
 239       fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord);
 240
 241       caloClusterPtr->fModule = recPointPtr->fModule;
 242       caloClusterPtr->fGlobalPos[0] =  globalCoord.fX;
 243       caloClusterPtr->fGlobalPos[1] =  globalCoord.fY;
 244       caloClusterPtr->fGlobalPos[2] =  globalCoord.fZ;
 245
 246       HLTDebug("Cluster local position: x = %f, z = %f, module = %d", recPointPtr->fX, recPointPtr->fZ, recPointPtr->fModule);
 247       HLTDebug("Cluster global position: x = %f, y = %f, z = %f", globalCoord.fX, globalCoord.fY, globalCoord.fZ);
 248
 249       //caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
 250       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
 251
 252       caloClusterPtr->fClusterType = fClusterType;
 253 //      Int_t tmpSize = 0;//totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 254
 255       //TODO remove hardcoded 10;
 256       memset(caloClusterPtr->fTracksMatched, 0xff, sizeof(Int_t)*10);
 257
 258       //Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 259       UInt_t tmpSize = (caloClusterPtr->fNCells-1)*sizeof(AliHLTCaloCellDataStruct);
 260
 261       if((availableSize - totSize)  < tmpSize)
 262       {
 263          HLTError("Out of buffer, available size is: %d, total size used: %d, extra size needed: %d", availableSize, totSize, tmpSize);
 264          return -ENOBUFS;
 265       }
 266
 267       Int_t *digitIndexPtr = &(recPointPtr->fDigits);
 268       Int_t id = 0;
 269
 270        AliHLTCaloCellDataStruct *cellPtr = &(caloClusterPtr->fCaloCells);
 271
 272       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
 273         {
 274            digitPtr = &(fDigitDataArray[*digitIndexPtr]);
 275            fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
 276
 277           cellPtr->fCellsAbsId= id;
 278           cellPtr->fCellsAmpFraction = digitPtr->fEnergy/recPointPtr->fAmp;
 279           //printf("Cell ID pointer: %x\n", cellIDPtr);
 280          //printf("Cell Amp Pointer: %x\n", cellAmpFracPtr);
 281          //printf("Cell pos: x = %d, z = %d\n", digitPtr->fX, digitPtr->fZ);
 282 //       printf("Cell ID: %d\n", cellPtr->fCellsAbsId);
 283          //printf("Cell Amp: %f, pointer: %x\n", *cellAmpFracPtr, cellAmpFracPtr);
 284           cellPtr++;
 285           digitIndexPtr++;
 286
 287         }
 288
 289       totSize += tmpSize;
 290
 291       if(fRecoParamsPtr)
 292       {
 293          caloClusterPtr->fEnergy = fRecoParamsPtr->GetCorrectedEnergy(recPointPtr->fAmp);
 294       }
 295       else
 296       {
 297          caloClusterPtr->fEnergy = recPointPtr->fAmp;
 298       }
 299
 300       HLTDebug("Cluster global position: x = %f, y = %f, z = %f, energy: %f, number of cells: %d, cluster pointer: %x", globalCoord.fX, globalCoord.fY, globalCoord.fZ, caloClusterPtr->fEnergy, caloClusterPtr->fNCells,  caloClusterPtr);
 301
 302       if(fDoClusterFit)
 303         {
 304           FitCluster(recPointPtr);
 305         }
 306       else
 307         {
 308           caloClusterPtr->fDispersion = 0;
 309           caloClusterPtr->fFitQuality = 0;
 310           caloClusterPtr->fM20 = 0;
 311           caloClusterPtr->fM02 = 0;
 312
 313         }
 314       if(fHaveCPVInfo)
 315         {
 316           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
 317         }
 318       else
 319         {
 320           caloClusterPtr->fEmcCpvDistance = -1;
 321         }
 322       if(fDoPID)
 323         {
 324           DoParticleIdentification(caloClusterPtr);
 325         }
 326       else
 327         {
 328           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
 329             {
 330               caloClusterPtr->fPID[k] = 0;
 331             }
 332         }
 333       if(fHaveDistanceToBadChannel)
 334         {
 335           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
 336         }
 337       else
 338         {
 339           caloClusterPtr->fDistanceToBadChannel = -1;
 340         }
 341
 342       caloClusterPtr->fClusterType = fClusterType;
 343           //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
 344       //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 345       //printf("CaloClusterPtr: %x, energy: %f\n", caloClusterPtr, caloClusterPtr->fEnergy);
 346
 347       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
 348       //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
 349       fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
 350       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
 351       //digitPtr = &(recPointPtr->fDigits);
 352     }
 353
 354 return fNRecPoints;
 355
 356 }
 357