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