]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/CALO/AliHLTCaloClusterAnalyser.cxx
RWGCF converted to native cmake
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterAnalyser.cxx
index a2ac118da908901ab1ba22d99b9e5435d9f2cd13..d38841bbab335f624c8969ac77724839fc461b71 100644 (file)
@@ -39,7 +39,7 @@
 #include "TVector3.h"
 #include "TH1F.h"
 #include "TFile.h"
-#include "AliHLTCaloClusterizer.h"
+#include "AliHLTCaloRecoParamHandler.h"
 
 ClassImp(AliHLTCaloClusterAnalyser);
 
@@ -57,7 +57,10 @@ AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
   fDoPID(false),
   fHaveDistanceToBadChannel(false),
   fGeometry(0),
-  fClusterType(AliESDCaloCluster::kPHOSCluster)
+  fClusterType(AliVCluster::kPHOSNeutral),
+  fRecoParamsPtr(0),
+  fCutOnSingleCellClusters(false),
+  fSingleCellEnergyCut(0.5)
 {
   //See header file for documentation
 }
@@ -82,7 +85,7 @@ AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPo
 }
 
 void 
-AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits) 
+AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct **digits) 
 { 
 //   AliHLTCaloClusterizer cl("PHOS");
   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
@@ -110,34 +113,55 @@ AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
       //      digit = &(recPoint->fDigits);
 
       Int_t *digitIndexPtr = &(recPoint->fDigits);
-
-      for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
-       {
-
-         digit = &(fDigitDataArray[*digitIndexPtr]);
-
-         xi = digit->fX+0.5;
-         zi = digit->fZ+0.5;
-
-         if (recPoint->fAmp > 0 && digit->fEnergy > 0) 
-           {
-             Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
-             x    += xi * w ;
-             z    += zi * w ;
-             wtot += w ;
-           }
-         digitIndexPtr++;
-       }
-
-      if (wtot>0) 
-       {
-         recPoint->fX = x/wtot ;
-         recPoint->fZ = z/wtot ;
-       }
-      else
-       {
-         recPoint->fAmp = 0;
-       }
+       wtot = 0;
+       x = 0; 
+       z = 0;
+
+/*       Float_t maxAmp = 0;
+       Int_t maxX = 0;
+       Int_t maxZ = 0;*/
+       if (fDigitDataArray[*digitIndexPtr])
+
+        for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
+          {
+            
+            digit = fDigitDataArray[*digitIndexPtr];
+            
+            xi = digit->fX;
+            zi = digit->fZ;
+            
+            //xi = digit->fX+0.5;
+            //zi = digit->fZ+0.5;
+            
+            if (recPoint->fAmp > 0 && digit->fEnergy > 0) 
+              {
+                Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
+                x    += xi * w ;
+                z    += zi * w ;
+                wtot += w ;
+                /*           if(digit->fEnergy > maxAmp)
+                             {
+                             maxAmp = digit->fEnergy;
+                             maxX = digit->fX;// + 0.5;
+                             maxZ = digit->fZ;// + 0.5;
+                             }*/
+              }
+            digitIndexPtr++;
+          }
+       
+       if (wtot>0) 
+        {
+          recPoint->fX = x/wtot ;
+          recPoint->fZ = z/wtot ;
+        }
+       else
+        {
+          recPoint->fX = -9999;
+          recPoint->fZ =-9999;
+          // no good crashes depth with FP exception
+          //recPoint->fAmp = 0;
+        }
+//     printf("Max digit: E = %f, x = %d, z= %d, cluster: E = %f, x = %f, z = %f\n" , maxAmp, maxX, maxZ, recPoint->fAmp, recPoint->fX, recPoint->fZ);
     }
   return 0;
 }
@@ -170,13 +194,17 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
 {
   //See header file for documentation
 
+   if(fRecoParamsPtr)
+   {
+      fLogWeight = fRecoParamsPtr->GetLogWeight();
+   }
    
-  totSize += sizeof(AliHLTCaloClusterDataStruct);
   fNRecPoints = nRecPoints;
 
   if(fGeometry == 0)
   {
      HLTError("No geometry object is initialised, creation of clusters stopped");
+     return  -1;
   }
 
   CalculateCenterOfGravity();
@@ -185,34 +213,55 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
   AliHLTCaloDigitDataStruct* digitPtr = 0;
 
   AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
-  UShort_t* cellIDPtr = 0;
-  Float_t* cellAmpFracPtr = 0;;
-  
-//  Int_t id = -1;
+
+  //  Int_t id = -1;
   TVector3 globalPos;
 
   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
     {
       if((availableSize - totSize)  < sizeof(AliHLTCaloClusterDataStruct))
       {
-        HLTError("Out of buffer");
+        HLTError("Out of buffer: available size is: %d, total size used: %d", availableSize, totSize);
         return -ENOBUFS;
       }
       
-      caloClusterPtr = fCaloClusterDataPtr;
-     
-      cellIDPtr = &(caloClusterPtr->fCellsAbsId);
-      cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
-     
       AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
+
+      if(fCutOnSingleCellClusters && recPointPtr->fAmp > fSingleCellEnergyCut && recPointPtr->fMultiplicity == 1) continue;
+      
+      totSize += sizeof(AliHLTCaloClusterDataStruct);
+      
+      caloClusterPtr = fCaloClusterDataPtr;
+      caloClusterPtr->fChi2 = 0;
+      caloClusterPtr->fClusterType = kUndef;
+      caloClusterPtr->fDispersion = 0;
+      caloClusterPtr->fDistanceToBadChannel = 0;
+      caloClusterPtr->fDistToBadChannel = 0;
+      caloClusterPtr->fEmcCpvDistance = 0;
+      caloClusterPtr->fEnergy = 0;
+      caloClusterPtr->fFitQuality = 0;
+      caloClusterPtr->fID = 0;
+      caloClusterPtr->fM02 = 0;
+      caloClusterPtr->fM20 = 0;
+      caloClusterPtr->fNCells = 0;
+      caloClusterPtr->fNExMax = 0;
+      caloClusterPtr->fTOF = 0;
+      caloClusterPtr->fTrackDx = -999;
+      caloClusterPtr->fTrackDz = -999;
       
       AliHLTCaloGlobalCoordinate globalCoord;
-      fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord);
 
+      // 0 = assume photon
+      fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord, 0);
+
+      caloClusterPtr->fModule = recPointPtr->fModule;
       caloClusterPtr->fGlobalPos[0] =  globalCoord.fX;
       caloClusterPtr->fGlobalPos[1] =  globalCoord.fY;
       caloClusterPtr->fGlobalPos[2] =  globalCoord.fZ;
 
+      HLTDebug("Cluster local position: x = %f, z = %f, module = %d", recPointPtr->fX, recPointPtr->fZ, recPointPtr->fModule);
+      HLTDebug("Cluster global position: x = %f, y = %f, z = %f", globalCoord.fX, globalCoord.fY, globalCoord.fZ);
+      
       //caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
 
@@ -222,38 +271,54 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
       //TODO remove hardcoded 10; 
       memset(caloClusterPtr->fTracksMatched, 0xff, sizeof(Int_t)*10);
 
-      Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
+      //Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
+      UInt_t tmpSize = (caloClusterPtr->fNCells-1)*sizeof(AliHLTCaloCellDataStruct);
 
       if((availableSize - totSize)  < tmpSize)
       {
-        HLTError("Out of buffer");
+        HLTError("Out of buffer, available size is: %d, total size used: %d, extra size needed: %d", availableSize, totSize, tmpSize);
         return -ENOBUFS;
       }
+      
       Int_t *digitIndexPtr = &(recPointPtr->fDigits);
       Int_t id = 0;
 
+       AliHLTCaloCellDataStruct *cellPtr = &(caloClusterPtr->fCaloCells);
+      Float_t maxTime = 0; //time of maximum amplitude cell is assigned to cluster 
       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
        {
-          digitPtr = &(fDigitDataArray[*digitIndexPtr]);
-          id++;
-          // //          fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
-         *cellIDPtr = id;
-         *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
+          digitPtr = fDigitDataArray[*digitIndexPtr];
+          fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
+          
+         cellPtr->fCellsAbsId= id;
+         cellPtr->fCellsAmpFraction = digitPtr->fEnergy/recPointPtr->fAmp;
+               if(digitPtr->fTime > maxTime) maxTime = digitPtr->fTime; 
          //printf("Cell ID pointer: %x\n", cellIDPtr);
         //printf("Cell Amp Pointer: %x\n", cellAmpFracPtr);
-        //printf("Digit ID: %d\n", *cellIDPtr);
-         digitPtr++;
-         //cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t)); 
-         cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Float_t) + sizeof(Short_t)); 
-         cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Short_t) + sizeof(Float_t));
-         
+        //printf("Cell pos: x = %d, z = %d\n", digitPtr->fX, digitPtr->fZ);
+//      printf("Cell ID: %d\n", cellPtr->fCellsAbsId);
+        //printf("Cell Amp: %f, pointer: %x\n", *cellAmpFracPtr, cellAmpFracPtr);
+         cellPtr++;
          digitIndexPtr++;
          
        }
 
       totSize += tmpSize;
 
-      caloClusterPtr->fEnergy = recPointPtr->fAmp;
+      if(fRecoParamsPtr)
+      {
+        caloClusterPtr->fEnergy = fRecoParamsPtr->GetCorrectedEnergy(recPointPtr->fAmp);
+      }
+      else
+      {
+        caloClusterPtr->fEnergy = recPointPtr->fAmp;
+      }
+      
+      // Set the time of the maximum digit as cluster time
+      //caloClusterPtr->fTOF = recPointPtr->fTime;
+                       caloClusterPtr->fTOF = maxTime;
+      
+      HLTDebug("Cluster global position: x = %f, y = %f, z = %f, energy: %f, time: %f, number of cells: %d, cluster pointer: %x", globalCoord.fX, globalCoord.fY, globalCoord.fZ, caloClusterPtr->fEnergy, caloClusterPtr->fTOF, caloClusterPtr->fNCells,  caloClusterPtr);
 
       if(fDoClusterFit)
        {
@@ -281,7 +346,7 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
        }
       else
        {
-         for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
+         for(Int_t k = 0; k < AliPID::kSPECIESCN; k++)
            {
              caloClusterPtr->fPID[k] = 0;
            }
@@ -302,7 +367,7 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
       
       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
       //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
-      fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
+      fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
       //digitPtr = &(recPointPtr->fDigits);  
     }