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