]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliCaloRawAnalyzerPeakFinder.cxx
EbE q vector correlations
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerPeakFinder.cxx
index f839260957e58f5c14776fecdf90af3d16af2bd4..e8bc413026a5babdcd1e56af0b941c1dce375363 100644 (file)
@@ -1,11 +1,12 @@
+// -*- mode: c++ -*-
 /**************************************************************************
- * This file is property of and copyright by the Experimental Nuclear     *
- * Physics Group, Dep. of Physics                                         *
- * University of Oslo, Norway, 2007                                       *
+ * This file is property of and copyright by                              *
+ * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
+ *                                                                        *
+ * Primary Author: Per Thomas Hille  <perthomas.hille@yale.edu>           *
  *                                                                        *
- * Author: Per Thomas Hille <perthi@fys.uio.no> for the ALICE HLT Project.*
  * Contributors are mentioned in the code where appropriate.              *
- * Please report bugs to perthi@fys.uio.no                                *
+ * Please report bugs to   perthomas.hille@yale.edu                       *
  *                                                                        *
  * Permission to use, copy, modify and distribute this software and its   *
  * documentation strictly for non-commercial purposes is hereby granted   *
 // The amplitude is extracted  as a
 // weighted sum of the samples using the 
 // best possible weights.
-// The wights is calculated only once and the
-// Actual extraction of amplitude and peak position
-// Is done with a simple vector multiplication, allowing for
-// Extreemely fast computations. 
+// The weights are calculated only once and the
+// actual extraction of amplitude and peak position
+// is done with a simple vector multiplication, allowing
+// extremely fast computations.
 
 #include "AliCaloRawAnalyzerPeakFinder.h"
 #include "AliCaloBunchInfo.h"
 #include "AliCaloFitResults.h"
-#include <iostream>
-#include "unistd.h"
 #include "TMath.h"
 #include "AliLog.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "TFile.h"
+#include "AliCaloPeakFinderVectors.h"
+#include <iostream>
 
 using namespace std;
 
-AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer(), 
-                                                               fTof(0), 
-                                                               fAmp(0)
-{
-  //comment
+ClassImp( AliCaloRawAnalyzerPeakFinder )
 
-  fNsampleCut = 5;
 
-  for(int i=0; i < MAXSTART; i++)
+AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),  
+                                                             fPeakFinderVectors(0),
+                                                             fRunOnAlien(false),
+                                                             fIsInitialized(false)
+{
+  // Ctor
+  fAlgo= Algo::kPeakFinder;
+  fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
+  ResetVectors();
+  LoadVectorsOCDB();
+}
+
+void  
+AliCaloRawAnalyzerPeakFinder::ResetVectors()
+{
+  //As name implies
+  for(int i=0; i < PF::MAXSTART; i++)
     {
-      for(int j=0; j < SAMPLERANGE; j++ )
+      for(int j=0; j < PF::SAMPLERANGE; j++ )
        {
-         fPFAmpVectors[i][j] = new double[100];
-         fPFTofVectors[i][j] = new double[100];
-         
          for(int k=0; k < 100; k++ )
            {
              fPFAmpVectors[i][j][k] = 0; 
              fPFTofVectors[i][j][k] = 0;
+             fPFAmpVectorsCoarse[i][j][k] = 0;
+             fPFTofVectorsCoarse[i][j][k] = 0; 
            }
        }
     }
-  LoadVectors();
 }
 
 
-AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
+Double_t  
+AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
 {
-  //comment
-  for(int i=0; i < MAXSTART; i++)
+  // Fisrt (coarce) estimate of Amplitude using the Peak-Finder.
+  // The output of the first iteration is sued to select vectors 
+  // for the second iteration.
+
+  Double_t tmpTof = 0;
+  Double_t tmpAmp= 0;
+
+  for(int i=0; i < length; i++)
     {
-      for(int j=0; j < SAMPLERANGE; j++ )
-       {
-         delete[] fPFAmpVectors[i][j];
-         delete[] fPFTofVectors[i][j];
-       }
+      tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i]; 
+      tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i]; 
     }
+  
+  tmpTof = tmpTof / tmpAmp ;
+  return tmpTof;
 }
 
 
 AliCaloFitResults 
 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
 {
+  // Evaluation of amplitude and TOF
+  if( fIsInitialized == false )
+    {
+      cout << __FILE__ << ":" << __LINE__ << "ERROR, peakfinder vectors not loaded" << endl;
+      return  AliCaloFitResults(kInvalid, kInvalid);
+    }
+
   // Extracting the amplitude using the Peak-Finder algorithm
   // The amplitude is a weighted sum of the samples using 
   // optimum weights.
@@ -86,87 +113,163 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
   short maxamp; //Maximum amplitude
   fAmp = 0;
   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
+  
   if( index >= 0)
     {
       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
-      
-      if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )         // (maxamp - ped) > 900 = Close to saturation (use low gain then)
+      short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1); 
+      if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
        {
-         //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
-         return  AliCaloFitResults( maxamp, ped,  -1, maxf,   maxampindex, -1, -1 );
-       }
-      
-      int first;
-      int last;
-      
-      if ( maxf > fAmpCut )
+         return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
+       }            
+      else if ( maxf >= fAmpCut )
        {
-         SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex -  bunchvector.at(index).GetStartBin(), &first, &last);
+         int first = 0;
+         int last = 0;
+         short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();     
+         SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
          int nsamples =  last - first;
-         if( ( nsamples  )  >= fNsampleCut )
+
+         if( ( nsamples  )  >= fNsampleCut ) // no if statement needed really; keep for readability
            {
              int startbin = bunchvector.at(index).GetStartBin();  
-             int n = last -first;  
+             int n = last - first;  
              int pfindex = n - fNsampleCut; 
-             pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
+             pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
+             int dt =  maxampindex - startbin -2; 
+             int tmpindex = 0;
+             Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
+             
+             if( tmptof < -1 )
+               {
+                 tmpindex = 0;
+               }
+             else
+               if( tmptof  > -1 && tmptof < 100 )
+                 {
+                   tmpindex = 1;
+                 }
+               else
+                 {
+                   tmpindex = 2;
+                 }
+
+             double tof = 0;
+             for(int k=0; k < PF::SAMPLERANGE; k++   )
+               {
+                 tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];   
+               }
+             for( int i=0; i < PF::SAMPLERANGE; i++ )
+               {
+                 {
+                   
+                   fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i  +tmpindex -1 ];
+                 }
+               }
 
-             for(int i=0; i < SAMPLERANGE; i++ )
+             if( TMath::Abs(  (maxf - fAmp  )/maxf )  >   0.1 )
                {
-                 int dt =  maxampindex - startbin -2;
-                 double tmp =  fReversed[ dt  +i];  
-                 fAmp += fPFAmpVectors[0][pfindex][i]*tmp; 
+                 fAmp = maxf;
                }
              
-             return AliCaloFitResults( maxamp, ped , -1, fAmp, -1, -1, -1 ); 
+             tof = timebinOffset - 0.01*tof/fAmp - fL1Phase/TIMEBINWITH; // clock
+             Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
+             Int_t ndf = last - first - 1; // nsamples - 2
+             return AliCaloFitResults( maxamp, ped , Ret::kFitPar, fAmp, tof, 
+                                       timebinOffset, chi2, ndf,
+                                       Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );  
            }
          else
            {
-             return AliCaloFitResults( maxamp, ped , -1, maxf, -1, -1, -1 ); 
+             Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
+             Int_t ndf = last - first - 1; // nsamples - 2
+             return AliCaloFitResults( maxamp, ped , Ret::kCrude, maxf, timebinOffset,
+                                       timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
            }
-       }
+       } // ampcut
     }
-  //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
-
-  return  AliCaloFitResults(-1, -1);
+  return  AliCaloFitResults(kInvalid, kInvalid);
 }
 
 
-void 
-AliCaloRawAnalyzerPeakFinder::LoadVectors()
+void   
+AliCaloRawAnalyzerPeakFinder::CopyVectors( const AliCaloPeakFinderVectors *const pfv )
 {
-  //Read in the Peak finder vecors from file
-  for(int i = 0;  i < MAXSTART ; i++)
+  // As name implies
+  if ( pfv != 0)
     {
-      for( int j=0; j < SAMPLERANGE; j++)
+      for(int i = 0;  i < PF::MAXSTART ; i++)
        {
-         char filename[256];
-         int n = j+fNsampleCut;
-         double start = (double)i+0.5;
-         sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
-         FILE *fp = fopen(filename, "r");
-         
-         if(fp == 0)
+         for( int j=0; j < PF::SAMPLERANGE; j++)  
            {
-             AliFatal( Form( "could not open file: %s", filename ) );
+             pfv->GetVector( i, j, fPFAmpVectors[i][j] ,  fPFTofVectors[i][j],    
+                             fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j]  ); 
+
+             fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
+                                            fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
            }
-         else
+       }
+    }
+  else
+    {
+      AliFatal( "pfv = ZERO !!!!!!!");
+    } 
+}
+
+
+void   
+AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
+{
+  //Loading of Peak-Finder  vectors from the 
+  //Offline Condition Database  (OCDB)
+  AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
+  
+  if( entry != 0 )
+  {
+    //cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
+    //entry->PrintMetaData();
+    AliCaloPeakFinderVectors  *pfv = (AliCaloPeakFinderVectors *)entry->GetObject(); 
+    if( pfv == 0 )
+    {
+      cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
+    }
+    CopyVectors( pfv );
+    
+    if( pfv != 0 )
+    {
+      fIsInitialized = true;
+    }
+  }
+}
+
+
+void   
+AliCaloRawAnalyzerPeakFinder::WriteRootFile() const
+{ // Utility function to write Peak-Finder vectors to an root file
+  // The output is used to create an OCDB entry.
+  fPeakFinderVectors->PrintVectors();
+  TFile *f = new TFile("peakfindervectors2.root",  "recreate" );
+  fPeakFinderVectors->Write();
+  f->Close();
+  delete f;
+}
+
+
+void 
+AliCaloRawAnalyzerPeakFinder::PrintVectors()
+{ // Utility function to write Peak-Finder vectors 
+  for(int i=0; i < 20; i++)
+    {
+      for( int j = 0; j < PF::MAXSTART; j ++ )
+       {
+         for( int k=0; k < PF::SAMPLERANGE; k++ )
            {
-             for(int m = 0; m < n ; m++ )
-               {
-                 fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
-               }
-             fscanf(fp, "\n");
-             
-             for(int m = 0; m < n ; m++ )
-               {
-                 fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
-               }
-             
-             fclose (fp);
+             cout << fPFAmpVectors[j][k][i] << "\t" ;
            }
        }
+      cout << endl;
     }
+  cout << __FILE__ << ":" << __LINE__ << ":.... DONE !!" << endl;
 }