]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/CALO/AliHLTCaloClusterAnalyser.cxx
Updates EvtGen Code
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterAnalyser.cxx
index 7467f23e3639c47c1bd890c73f2e574e6ea170a7..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,18 +57,17 @@ 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
-  fHist = new TH1F("cluster_energies", "cluster_energies", 200, 0, 20);
-  
 }
 
 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser() 
 {
-   TFile *file = TFile::Open("debug.root", "RECREATE");
-   fHist->Write();
-   file->Close();
+  // See header file for class documentation
 }
 
 void 
@@ -86,7 +85,7 @@ AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPo
 }
 
 void 
-AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits) 
+AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct **digits) 
 { 
 //   AliHLTCaloClusterizer cl("PHOS");
   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
@@ -105,11 +104,6 @@ AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
   Float_t zi = 0.;
 
   AliHLTCaloDigitDataStruct *digit = 0;
-  //AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
-   
-  
-  //AliHLTCaloClusterizer      cl("PHOS");
-  //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
 
   UInt_t iDigit = 0;
 
@@ -117,38 +111,57 @@ AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
     {
       AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
       //      digit = &(recPoint->fDigits);
-      cout << "COG: Rec point multiplicity: " << recPoint->fMultiplicity << ", rec point energy: " << recPoint->fAmp << endl;
-      Int_t *digitIndexPtr = &(recPoint->fDigits);
 
-      for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
-       {
-         //cout << "1. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
-
-         digit = &(fDigitDataArray[*digitIndexPtr]);
-         cout << "COG: 2. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", index pointer: " << digitIndexPtr<<  endl;
-         cout << ", dig Energy: " << digit->fEnergy << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;        
-         xi = digit->fX;
-         zi = digit->fZ;
-         //  cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
-         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++;
-       }
-      //cout << endl;
-      if (wtot>0) 
-       {
-         recPoint->fX = x/wtot ;
-         recPoint->fZ = z/wtot ;
-       }
-      else
-       {
-         recPoint->fAmp = 0;
-       }
+      Int_t *digitIndexPtr = &(recPoint->fDigits);
+       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;
 }
@@ -156,7 +169,7 @@ AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
 
 Int_t 
 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
-{
+{      
   //See header file for documentation
   return 0;
 }
@@ -181,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();
@@ -196,67 +213,113 @@ 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];
-      cout << "CA: rec point energy: " << recPointPtr->fAmp << ", rec point multiplicity: " << recPointPtr->fMultiplicity << endl;
-      //      cout << "Local Position (x:z:module): " << recPointPtr->fX << " : "<< recPointPtr->fZ << " : " << recPointPtr->fModule << endl;
+
+      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);
-      // cout << "Global Position (x:y:z): " << globalPos[0] << " : "<< globalPos[1] << " : " << globalPos[2] << endl << endl;
 
-      caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
-      caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
-      caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
+      // 0 = assume photon
+      fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord, 0);
 
-      caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
-      
-      Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
+      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;
+
+      caloClusterPtr->fClusterType = fClusterType;
+//      Int_t tmpSize = 0;//totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
+
+      //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));
+      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++)
        {
-/*       fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
-         *cellIDPtr = id;
-         *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
-         digitPtr++;
-         cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t)); 
-         cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t)); */
+          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("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;
-
-      fHist->Fill(caloClusterPtr->fEnergy);
-
-      cout << "CA: cluster energy: " << caloClusterPtr->fEnergy << endl;
-      cout << "CA: recpoint energy: " << recPointPtr->fAmp << endl;
+      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)
        {
          FitCluster(recPointPtr);
@@ -283,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;
            }
@@ -298,18 +361,18 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
        }
 
       caloClusterPtr->fClusterType = fClusterType;
-      //      totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);   
+         //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);   
       //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));   
-
+      //printf("CaloClusterPtr: %x, energy: %f\n", caloClusterPtr, caloClusterPtr->fEnergy);
+      
       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
-      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
-
+      //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
+      fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
       //digitPtr = &(recPointPtr->fDigits);  
     }
-  //  cout << "CA: Energy End: " << fCaloClusterDataPtr->fEnergy << endl;
-  //cout << "CA totSize: " << totSize << endl;
-  return fNRecPoints;
+
+return fNRecPoints;
 
 }