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 #ifndef HAVE_NOT_ALIVCLUSTER // backward compatibility for r42844
  61   fClusterType(AliVCluster::kPHOSNeutral),
  62 #else
  63   fClusterType(AliESDCaloCluster::kPHOSCluster),
  64 #endif
  65   fRecoParamsPtr(0),
  66   fCutOnSingleCellClusters(false),
  67   fSingleCellEnergyCut(0.5)
  68 {
  69   //See header file for documentation
  70 }
  71
  72 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
  73 {
  74   // See header file for class documentation
  75 }
  76
  77 void
  78 AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
  79 {
  80   //see header file for documentation
  81   fCaloClusterDataPtr = caloClusterDataPtr;
  82 }
  83
  84 void
  85 AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
  86 {
  87   fRecPointArray = recPointDataPtr;
  88   fNRecPoints = nRecPoints;
  89 }
  90
  91 void
  92 AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct **digits)
  93 {
  94 //   AliHLTCaloClusterizer cl("PHOS");
  95   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
  96    fDigitDataArray = digits;
  97    //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
  98 }
  99
 100 Int_t
 101 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
 102 {
 103   //see header file for documentation
 104   Float_t wtot = 0.;
 105   Float_t x = 0.;
 106   Float_t z = 0.;
 107   Float_t xi = 0.;
 108   Float_t zi = 0.;
 109
 110   AliHLTCaloDigitDataStruct *digit = 0;
 111
 112   UInt_t iDigit = 0;
 113
 114   for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
 115     {
 116       AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
 117       //      digit = &(recPoint->fDigits);
 118
 119       Int_t *digitIndexPtr = &(recPoint->fDigits);
 120        wtot = 0;
 121        x = 0;
 122        z = 0;
 123
 124 /*       Float_t maxAmp = 0;
 125        Int_t maxX = 0;
 126        Int_t maxZ = 0;*/
 127        if (fDigitDataArray[*digitIndexPtr])
 128
 129          for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 130            {
 131
 132              digit = fDigitDataArray[*digitIndexPtr];
 133
 134              xi = digit->fX;
 135              zi = digit->fZ;
 136
 137              //xi = digit->fX+0.5;
 138              //zi = digit->fZ+0.5;
 139
 140              if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 141                {
 142                  Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 143                  x    += xi * w ;
 144                  z    += zi * w ;
 145                  wtot += w ;
 146                  /*           if(digit->fEnergy > maxAmp)
 147                               {
 148                               maxAmp = digit->fEnergy;
 149                               maxX = digit->fX;// + 0.5;
 150                               maxZ = digit->fZ;// + 0.5;
 151                               }*/
 152                }
 153              digitIndexPtr++;
 154            }
 155
 156        if (wtot>0)
 157          {
 158            recPoint->fX = x/wtot ;
 159            recPoint->fZ = z/wtot ;
 160          }
 161        else
 162          {
 163            recPoint->fX = -9999;
 164            recPoint->fZ =-9999;
 165            // no good crashes depth with FP exception
 166            //recPoint->fAmp = 0;
 167          }
 168 //      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);
 169     }
 170   return 0;
 171 }
 172
 173
 174 Int_t
 175 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
 176 {
 177   //See header file for documentation
 178   return 0;
 179 }
 180
 181 Int_t
 182 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
 183 {
 184   //See header file for documentation
 185   return 0;
 186 }
 187
 188
 189 Int_t
 190 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
 191 {
 192   //See header file for documentation
 193   return 0;
 194 }
 195
 196 Int_t
 197 AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
 198 {
 199   //See header file for documentation
 200
 201    if(fRecoParamsPtr)
 202    {
 203       fLogWeight = fRecoParamsPtr->GetLogWeight();
 204    }
 205
 206   fNRecPoints = nRecPoints;
 207
 208   if(fGeometry == 0)
 209   {
 210      HLTError("No geometry object is initialised, creation of clusters stopped");
 211      return  -1;
 212   }
 213
 214   CalculateCenterOfGravity();
 215
 216   //  AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
 217   AliHLTCaloDigitDataStruct* digitPtr = 0;
 218
 219   AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
 220
 221   //  Int_t id = -1;
 222   TVector3 globalPos;
 223
 224   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
 225     {
 226       if((availableSize - totSize)  < sizeof(AliHLTCaloClusterDataStruct))
 227       {
 228          HLTError("Out of buffer: available size is: %d, total size used: %d", availableSize, totSize);
 229          return -ENOBUFS;
 230       }
 231
 232       AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
 233
 234       if(fCutOnSingleCellClusters && recPointPtr->fAmp > fSingleCellEnergyCut && recPointPtr->fMultiplicity == 1) continue;
 235
 236       totSize += sizeof(AliHLTCaloClusterDataStruct);
 237
 238       caloClusterPtr = fCaloClusterDataPtr;
 239       caloClusterPtr->fChi2 = 0;
 240       caloClusterPtr->fClusterType = kUndef;
 241       caloClusterPtr->fDispersion = 0;
 242       caloClusterPtr->fDistanceToBadChannel = 0;
 243       caloClusterPtr->fDistToBadChannel = 0;
 244       caloClusterPtr->fEmcCpvDistance = 0;
 245       caloClusterPtr->fEnergy = 0;
 246       caloClusterPtr->fFitQuality = 0;
 247       caloClusterPtr->fID = 0;
 248       caloClusterPtr->fM02 = 0;
 249       caloClusterPtr->fM20 = 0;
 250       caloClusterPtr->fNCells = 0;
 251       caloClusterPtr->fNExMax = 0;
 252       caloClusterPtr->fTOF = 0;
 253       caloClusterPtr->fTrackDx = -999;
 254       caloClusterPtr->fTrackDz = -999;
 255
 256       AliHLTCaloGlobalCoordinate globalCoord;
 257
 258       // 0 = assume photon
 259       fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord, 0);
 260
 261       caloClusterPtr->fModule = recPointPtr->fModule;
 262       caloClusterPtr->fGlobalPos[0] =  globalCoord.fX;
 263       caloClusterPtr->fGlobalPos[1] =  globalCoord.fY;
 264       caloClusterPtr->fGlobalPos[2] =  globalCoord.fZ;
 265
 266       HLTDebug("Cluster local position: x = %f, z = %f, module = %d", recPointPtr->fX, recPointPtr->fZ, recPointPtr->fModule);
 267       HLTDebug("Cluster global position: x = %f, y = %f, z = %f", globalCoord.fX, globalCoord.fY, globalCoord.fZ);
 268
 269       //caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
 270       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
 271
 272       caloClusterPtr->fClusterType = fClusterType;
 273 //      Int_t tmpSize = 0;//totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 274
 275       //TODO remove hardcoded 10;
 276       memset(caloClusterPtr->fTracksMatched, 0xff, sizeof(Int_t)*10);
 277
 278       //Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 279       UInt_t tmpSize = (caloClusterPtr->fNCells-1)*sizeof(AliHLTCaloCellDataStruct);
 280
 281       if((availableSize - totSize)  < tmpSize)
 282       {
 283          HLTError("Out of buffer, available size is: %d, total size used: %d, extra size needed: %d", availableSize, totSize, tmpSize);
 284          return -ENOBUFS;
 285       }
 286
 287       Int_t *digitIndexPtr = &(recPointPtr->fDigits);
 288       Int_t id = 0;
 289
 290        AliHLTCaloCellDataStruct *cellPtr = &(caloClusterPtr->fCaloCells);
 291       Float_t maxTime = 0; //time of maximum amplitude cell is assigned to cluster
 292       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
 293         {
 294            digitPtr = fDigitDataArray[*digitIndexPtr];
 295            fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
 296
 297           cellPtr->fCellsAbsId= id;
 298           cellPtr->fCellsAmpFraction = digitPtr->fEnergy/recPointPtr->fAmp;
 299                 if(digitPtr->fTime > maxTime) maxTime = digitPtr->fTime;
 300           //printf("Cell ID pointer: %x\n", cellIDPtr);
 301          //printf("Cell Amp Pointer: %x\n", cellAmpFracPtr);
 302          //printf("Cell pos: x = %d, z = %d\n", digitPtr->fX, digitPtr->fZ);
 303 //       printf("Cell ID: %d\n", cellPtr->fCellsAbsId);
 304          //printf("Cell Amp: %f, pointer: %x\n", *cellAmpFracPtr, cellAmpFracPtr);
 305           cellPtr++;
 306           digitIndexPtr++;
 307
 308         }
 309
 310       totSize += tmpSize;
 311
 312       if(fRecoParamsPtr)
 313       {
 314          caloClusterPtr->fEnergy = fRecoParamsPtr->GetCorrectedEnergy(recPointPtr->fAmp);
 315       }
 316       else
 317       {
 318          caloClusterPtr->fEnergy = recPointPtr->fAmp;
 319       }
 320
 321       // Set the time of the maximum digit as cluster time
 322       //caloClusterPtr->fTOF = recPointPtr->fTime;
 323                         caloClusterPtr->fTOF = maxTime;
 324
 325       HLTDebug("Cluster global position: x = %f, y = %f, z = %f, energy: %f, time: %f, number of cells: %d, cluster pointer: %x", globalCoord.fX, globalCoord.fY, globalCoord.fZ, caloClusterPtr->fEnergy, caloClusterPtr->fTOF, caloClusterPtr->fNCells,  caloClusterPtr);
 326
 327       if(fDoClusterFit)
 328         {
 329           FitCluster(recPointPtr);
 330         }
 331       else
 332         {
 333           caloClusterPtr->fDispersion = 0;
 334           caloClusterPtr->fFitQuality = 0;
 335           caloClusterPtr->fM20 = 0;
 336           caloClusterPtr->fM02 = 0;
 337
 338         }
 339       if(fHaveCPVInfo)
 340         {
 341           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
 342         }
 343       else
 344         {
 345           caloClusterPtr->fEmcCpvDistance = -1;
 346         }
 347       if(fDoPID)
 348         {
 349           DoParticleIdentification(caloClusterPtr);
 350         }
 351       else
 352         {
 353           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
 354             {
 355               caloClusterPtr->fPID[k] = 0;
 356             }
 357         }
 358       if(fHaveDistanceToBadChannel)
 359         {
 360           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
 361         }
 362       else
 363         {
 364           caloClusterPtr->fDistanceToBadChannel = -1;
 365         }
 366
 367       caloClusterPtr->fClusterType = fClusterType;
 368           //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
 369       //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 370       //printf("CaloClusterPtr: %x, energy: %f\n", caloClusterPtr, caloClusterPtr->fEnergy);
 371
 372       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
 373       //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
 374       fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
 375       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
 376       //digitPtr = &(recPointPtr->fDigits);
 377     }
 378
 379 return fNRecPoints;
 380
 381 }
 382