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