HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx

   1 // $Id$
   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   AliHLTPHOSClusterAnalyser.cxx
  20  * @author Oystein Djuvsland
  21  * @date
  22  * @brief  Cluster analyser for PHOS 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 "AliHLTPHOSClusterAnalyser.h"
  32 #include "AliHLTPHOSRecPointHeaderStruct.h"
  33 #include "AliHLTPHOSRecPointDataStruct.h"
  34 #include "AliHLTCaloClusterDataStruct.h"
  35 #include "AliHLTPHOSPhysicsAnalyzer.h"
  36 #include "AliPHOSGeoUtils.h"
  37 #include "AliESDCaloCluster.h"
  38 #include "TMath.h"
  39 #include "TVector3.h"
  40
  41 ClassImp(AliHLTPHOSClusterAnalyser);
  42
  43 AliHLTPHOSClusterAnalyser::AliHLTPHOSClusterAnalyser() :
  44   //  AliHLTPHOSBase(),
  45   fLogWeight(4.5),
  46   fRecPointDataPtr(0),
  47   fNRecPoints(0),
  48   fCaloClusterDataPtr(0),
  49   fCaloClusterHeaderPtr(0),
  50   fPHOSGeometry(0),
  51   fDoClusterFit(false),
  52   fHaveCPVInfo(false),
  53   fDoPID(false),
  54   fHaveDistanceToBadChannel(false)
  55 {
  56   //See header file for documentation
  57 }
  58
  59 AliHLTPHOSClusterAnalyser::~AliHLTPHOSClusterAnalyser()
  60 {
  61 }
  62
  63 void
  64 AliHLTPHOSClusterAnalyser::SetCaloClusterDataPtr(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
  65 {
  66   //see header file for documentation
  67   fCaloClusterDataPtr = caloClusterDataPtr;
  68 }
  69 void
  70 AliHLTPHOSClusterAnalyser::SetRecPointDataPtr(AliHLTPHOSRecPointHeaderStruct *recPointDataPtr, AliHLTPHOSDigitHeaderStruct *digitHeaderPtr)
  71 {
  72   fNRecPoints = recPointDataPtr->fNRecPoints;
  73   fRecPointDataPtr = reinterpret_cast<AliHLTPHOSRecPointDataStruct*>(reinterpret_cast<Char_t*>(recPointDataPtr)+sizeof(AliHLTPHOSRecPointHeaderStruct));
  74   fDigitHeaderPtr = digitHeaderPtr;
  75 }
  76
  77 Int_t
  78 AliHLTPHOSClusterAnalyser::CalculateCenterOfGravity()
  79 {
  80   //see header file for documentation
  81   Float_t wtot = 0.;
  82   Float_t x = 0.;
  83   Float_t z = 0.;
  84   Float_t xi = 0.;
  85   Float_t zi = 0.;
  86
  87   AliHLTPHOSDigitDataStruct *digit = 0;
  88   //AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
  89
  90   AliHLTPHOSRecPointDataStruct *recPoint = fRecPointDataPtr;
  91
  92   UInt_t iDigit = 0;
  93
  94   for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
  95     {
  96       digit = reinterpret_cast<AliHLTPHOSDigitDataStruct*>(reinterpret_cast<Int_t>(fDigitHeaderPtr) + recPoint->fStartDigitOffset);
  97       AliHLTPHOSDigitReader reader;
  98       reader->SetCurrentDigit(digit);
  99       while(digit)
 100         {
 101           xi = digit->fX;
 102           zi = digit->fZ;
 103           //  cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
 104           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 105             {
 106               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 107               x    += xi * w ;
 108               z    += zi * w ;
 109               wtot += w ;
 110             }
 111           digit = reader->NextDigit();
 112         }
 113       if (wtot>0)
 114         {
 115           recPoint->fX = x/wtot ;
 116           recPoint->fZ = z/wtot ;
 117         }
 118       else
 119         {
 120           recPoint->fAmp = 0;
 121         }
 122       recPoint = reinterpret_cast<AliHLTPHOSRecPointDataStruct*>(digit);
 123     }
 124   return 0;
 125 }
 126
 127
 128 Int_t
 129 AliHLTPHOSClusterAnalyser::CalculateRecPointMoments()
 130 {
 131   //See header file for documentation
 132   return 0;
 133 }
 134
 135 Int_t
 136 AliHLTPHOSClusterAnalyser::CalculateClusterMoments(AliHLTPHOSRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
 137 {
 138   //See header file for documentation
 139   return 0;
 140 }
 141
 142
 143 Int_t
 144 AliHLTPHOSClusterAnalyser::DeconvoluteClusters()
 145 {
 146   //See header file for documentation
 147   return 0;
 148 }
 149
 150 Int_t
 151 AliHLTPHOSClusterAnalyser::CreateClusters(UInt_t availableSize, UInt_t& totSize)
 152 {
 153   //See header file for documentation
 154
 155   UInt_t maxClusterSize = sizeof(AliHLTCaloClusterDataStruct) + (6 << 7); //Reasonable estimate... (6 = sizeof(Short_t) + sizeof(Float_t)
 156
 157   AliHLTPHOSRecPointDataStruct* recPointPtr = fRecPointDataPtr;
 158   AliHLTPHOSDigitDataStruct* digitPtr = reinterpret_cast<AliHLTPHOSDigitDataStruct*>(reinterpret_cast<Int_t>(fDigitHeaderPtr) + recPoint->fStartDigitOffset);
 159
 160   AliHLTPHOSReader reader;
 161   reader.SetCurrentDigit(digitPtr);
 162
 163   AliHLTCaloClusterDataStruct* caloClusterPtr = fCaloClusterDataPtr;
 164   UShort_t* cellIDPtr = &(caloClusterPtr->fCellsAbsId);
 165   Float_t* cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
 166
 167   Int_t id = -1;
 168   TVector3 globalPos;
 169
 170   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
 171     {
 172
 173       if(availableSize < (totSize + maxClusterSize))
 174         {
 175           return -1; //Might get out of buffer, exiting
 176         }
 177       fPHOSGeometry->Local2Global(recPointPtr->fModule, recPointPtr->fX, recPointPtr->fZ, globalPos);
 178
 179       caloClusterPtr->fGlobalPos[0] = globalPos[0];
 180       caloClusterPtr->fGlobalPos[1] = globalPos[1];
 181       caloClusterPtr->fGlobalPos[2] = globalPos[2];
 182
 183       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
 184
 185       cellIDPtr = &(caloClusterPtr->fCellsAbsId);
 186       cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
 187
 188       while(digitPtr)
 189         {
 190           fPHOSGeometry->RelPosToAbsId((Int_t)(recPointPtr->fModule), (double)(digitPtr->fX), (double)(digitPtr->fZ), id);
 191           *cellIDPtr = id;
 192           *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
 193           cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t));
 194           cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t));
 195           digitPtr = reader.NextDigit()
 196         }
 197
 198       caloClusterPtr->fEnergy = recPointPtr->fAmp;
 199
 200       if(fDoClusterFit)
 201         {
 202           FitCluster(recPointPtr);
 203         }
 204       else
 205         {
 206           caloClusterPtr->fDispersion = 0;
 207           caloClusterPtr->fFitQuality = 0;
 208           caloClusterPtr->fM20 = 0;
 209           caloClusterPtr->fM02 = 0;
 210
 211         }
 212       if(fHaveCPVInfo)
 213         {
 214           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
 215         }
 216       else
 217         {
 218           caloClusterPtr->fEmcCpvDistance = -1;
 219         }
 220       if(fDoPID)
 221         {
 222           DoParticleIdentification(caloClusterPtr);
 223         }
 224       else
 225         {
 226           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
 227             {
 228               caloClusterPtr->fPID[k] = 0;
 229             }
 230         }
 231       if(fHaveDistanceToBadChannel)
 232         {
 233           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
 234         }
 235       else
 236         {
 237           caloClusterPtr->fDistanceToBadChannel = -1;
 238         }
 239
 240       caloClusterPtr->fClusterType = (AliESDCaloCluster::kPHOSCluster);
 241
 242       totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 243
 244       caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
 245       recPointPtr = reinterpret_cast<AliHLTPHOSRecPointDataStruct*>(digitPtr);
 246       digitPtr = &(recPointPtr->fDigits);
 247     }
 248
 249   return fNRecPoints;
 250
 251 }
 252