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 "AliHLTCaloClusterizer.h"
  41
  42 ClassImp(AliHLTCaloClusterAnalyser);
  43
  44 AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
  45   //  AliHLTCaloBase(),
  46   fLogWeight(4.5),
  47   fRecPointArray(0),
  48   fDigitDataArray(0),
  49   fNRecPoints(0),
  50   fCaloClusterDataPtr(0),
  51   fCaloClusterHeaderPtr(0),
  52   //fAnalyzerPtr(0),
  53   fDoClusterFit(false),
  54   fHaveCPVInfo(false),
  55   fDoPID(false),
  56   fHaveDistanceToBadChannel(false),
  57   fGeometry(0)
  58 {
  59   //See header file for documentation
  60 }
  61
  62 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
  63 {
  64 }
  65
  66 void
  67 AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
  68 {
  69   //see header file for documentation
  70   fCaloClusterDataPtr = caloClusterDataPtr;
  71 }
  72
  73 void
  74 AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
  75 {
  76   fRecPointArray = recPointDataPtr;
  77   fNRecPoints = nRecPoints;
  78 }
  79
  80 void
  81 AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits)
  82 {
  83 //   AliHLTCaloClusterizer cl("PHOS");
  84   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
  85    fDigitDataArray = digits;
  86    //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
  87 }
  88
  89 Int_t
  90 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
  91 {
  92   //see header file for documentation
  93   Float_t wtot = 0.;
  94   Float_t x = 0.;
  95   Float_t z = 0.;
  96   Float_t xi = 0.;
  97   Float_t zi = 0.;
  98
  99   AliHLTCaloDigitDataStruct *digit = 0;
 100   //AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 101
 102
 103   //AliHLTCaloClusterizer       cl("PHOS");
 104   //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
 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       cout << "COG: Rec point multiplicity: " << recPoint->fMultiplicity << ", rec point energy: " << recPoint->fAmp << endl;
 113       Int_t *digitIndexPtr = &(recPoint->fDigits);
 114
 115       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
 116         {
 117           //cout << "1. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
 118
 119           digit = &(fDigitDataArray[*digitIndexPtr]);
 120           cout << "COG: 2. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", index pointer: " << digitIndexPtr<<  endl;
 121           cout << ", dig Energy: " << digit->fEnergy << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
 122           xi = digit->fX;
 123           zi = digit->fZ;
 124           //  cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
 125           if (recPoint->fAmp > 0 && digit->fEnergy > 0)
 126             {
 127               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
 128               x    += xi * w ;
 129               z    += zi * w ;
 130               wtot += w ;
 131             }
 132           digitIndexPtr++;
 133         }
 134       //cout << endl;
 135       if (wtot>0)
 136         {
 137           recPoint->fX = x/wtot ;
 138           recPoint->fZ = z/wtot ;
 139         }
 140       else
 141         {
 142           recPoint->fAmp = 0;
 143         }
 144     }
 145   return 0;
 146 }
 147
 148
 149 Int_t
 150 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
 151 {
 152   //See header file for documentation
 153   return 0;
 154 }
 155
 156 Int_t
 157 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
 158 {
 159   //See header file for documentation
 160   return 0;
 161 }
 162
 163
 164 Int_t
 165 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
 166 {
 167   //See header file for documentation
 168   return 0;
 169 }
 170
 171 Int_t
 172 AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
 173 {
 174   //See header file for documentation
 175
 176   fNRecPoints = nRecPoints;
 177
 178   if(fGeometry == 0)
 179   {
 180      HLTError("No geometry object is initialised, creation of clusters stopped");
 181   }
 182
 183   CalculateCenterOfGravity();
 184
 185   UInt_t maxClusterSize = sizeof(AliHLTCaloClusterDataStruct) + (6 << 7); //Reasonable estimate... (6 = sizeof(Short_t) + sizeof(Float_t)
 186
 187   //  AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
 188   AliHLTCaloDigitDataStruct* digitPtr = 0;
 189
 190   AliHLTCaloClusterDataStruct* caloClusterPtr = fCaloClusterDataPtr;
 191   UShort_t* cellIDPtr = &(caloClusterPtr->fCellsAbsId);
 192   Float_t* cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
 193
 194   Int_t id = -1;
 195   TVector3 globalPos;
 196
 197   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
 198     {
 199 //      if(availableSize
 200       AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
 201       cout << "CA: rec point energy: " << recPointPtr->fAmp << ", rec point multiplicity: " << recPointPtr->fMultiplicity << endl;
 202       if(availableSize < (totSize + maxClusterSize))
 203         {
 204           return -1; //Might get out of buffer, exiting
 205         }
 206       //      cout << "Local Position (x:z:module): " << recPointPtr->fX << " : "<< recPointPtr->fZ << " : " << recPointPtr->fModule << endl;
 207
 208       AliHLTCaloGlobalCoordinate globalCoord;
 209       fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord);
 210       // cout << "Global Position (x:y:z): " << globalPos[0] << " : "<< globalPos[1] << " : " << globalPos[2] << endl << endl;
 211
 212       caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
 213       caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
 214       caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
 215
 216       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
 217
 218       cellIDPtr = &(caloClusterPtr->fCellsAbsId);
 219       cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
 220
 221       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
 222         {
 223          // fPHOSGeometry->RelPosToAbsId((Int_t)(recPointPtr->fModule), (double)(digitPtr->fX), (double)(digitPtr->fZ), id);
 224 /*        *cellIDPtr = id;
 225           *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
 226           digitPtr++;
 227           cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t));
 228           cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t)); */
 229         }
 230
 231       caloClusterPtr->fEnergy = recPointPtr->fAmp;
 232       cout << "CA: cluster energy: " << caloClusterPtr->fEnergy << endl;
 233       cout << "CA: recpoint energy: " << recPointPtr->fAmp << endl;
 234       if(fDoClusterFit)
 235         {
 236           FitCluster(recPointPtr);
 237         }
 238       else
 239         {
 240           caloClusterPtr->fDispersion = 0;
 241           caloClusterPtr->fFitQuality = 0;
 242           caloClusterPtr->fM20 = 0;
 243           caloClusterPtr->fM02 = 0;
 244
 245         }
 246       if(fHaveCPVInfo)
 247         {
 248           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
 249         }
 250       else
 251         {
 252           caloClusterPtr->fEmcCpvDistance = -1;
 253         }
 254       if(fDoPID)
 255         {
 256           DoParticleIdentification(caloClusterPtr);
 257         }
 258       else
 259         {
 260           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
 261             {
 262               caloClusterPtr->fPID[k] = 0;
 263             }
 264         }
 265       if(fHaveDistanceToBadChannel)
 266         {
 267           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
 268         }
 269       else
 270         {
 271           caloClusterPtr->fDistanceToBadChannel = -1;
 272         }
 273
 274       caloClusterPtr->fClusterType = (AliESDCaloCluster::kPHOSCluster);
 275       //      totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
 276       totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
 277
 278
 279       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
 280       caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
 281
 282       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
 283       ///digitPtr = &(recPointPtr->fDigits);
 284     }
 285   //  cout << "CA: Energy End: " << fCaloClusterDataPtr->fEnergy << endl;
 286   //cout << "CA totSize: " << totSize << endl;
 287   return fNRecPoints;
 288
 289 }
 290