]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliCaloRawAnalyzerPeakFinder.cxx
Added method for filling TProfile
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerPeakFinder.cxx
index cd5cf733cd4233e72aed91928cebc44099ddbaf9..c246278d6af65a937e7bdc7c7164f0fca4b08c11 100644 (file)
@@ -1,3 +1,4 @@
+// -*- mode: c++ -*-
 /**************************************************************************
  * This file is property of and copyright by                              *
  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
@@ -5,7 +6,7 @@
  * Primary Author: Per Thomas Hille  <perthomas.hille@yale.edu>           *
  *                                                                        *
  * Contributors are mentioned in the code where appropriate.              *
- * Please report bugs to p.t.hille@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   *
 #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;
 
+
 ClassImp( AliCaloRawAnalyzerPeakFinder )
 
-AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF")
-//    fTof(0), 
-//                                                           fAmp(0)
+
+AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),  
+                                                             fAmp(0),
+                                                             fPeakFinderVectors(0),
+                                                             fRunOnAlien(false),
+                                                             fIsInitialized(false)
 {
-  //comment
+  //Comment
+  fAlgo= Algo::kPeakFinder;
+  InitOCDB(fRunOnAlien);
+  fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
+  ResetVectors();
+  LoadVectorsOCDB();
+}
 
-  fNsampleCut = 5;
 
-  for(int i=0; i < MAXSTART; i++)
+void 
+AliCaloRawAnalyzerPeakFinder::InitOCDB(bool alien) const
+{
+  // Setting the default OCDB pathe depending on wether we work locally or on the GRID.
+  AliCDBManager::Instance()->SetDefaultStorage(  alien == true ? "alien://$ALICE_ROOT/OCDB" : "local://$ALICE_ROOT/OCDB");
+  AliCDBManager::Instance()->SetRun(100);
+}
+
+
+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];
-         fPFAmpVectorsCoarse[i][j] = new double[100];
-         fPFTofVectorsCoarse[i][j] = new double[100];
-
          for(int k=0; k < 100; k++ )
            {
              fPFAmpVectors[i][j][k] = 0; 
@@ -63,31 +83,22 @@ AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer
            }
        }
     }
-
-  LoadVectors();
-
 }
 
 
 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
 {
   //comment
-  for(int i=0; i < MAXSTART; i++)
-    {
-      for(int j=0; j < SAMPLERANGE; j++ )
-       {
-         delete[] fPFAmpVectors[i][j];
-         delete[] fPFTofVectors[i][j];
-         delete[] fPFAmpVectorsCoarse[i][j];
-         delete[] fPFTofVectorsCoarse[i][j];
-       }
-    }
 }
 
 
 Double_t  
 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
 {
+  // 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;
 
@@ -105,84 +116,50 @@ AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int
 AliCaloFitResults 
 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
 {
+  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.
 
   short maxampindex; //index of maximum amplitude
   short maxamp; //Maximum amplitude
-  //  fAmp = 0;
-
-  
-  fAmpA[0] = 0;
-  fAmpA[1] = 0;
-  fAmpA[2] = 0;
-  
-  
-  //  cout << __FILE__ << __LINE__ << "\tendbin = " << bunchvector.at(index).GetEndBin()  <<  "\tstartbin = " << bunchvector.at(index).GetStartBin()  << endl;
-  
+  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);
          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 pfindex = n - fNsampleCut; 
-             pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
+             pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
 
-             short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1); 
-            
              int dt =  maxampindex - startbin -2; 
-
-             //    cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
-           
-        
-             //     Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
-             
-             //              cout << __FILE__ << __LINE__ << ",  dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
-
-             for( int i=0; i < SAMPLERANGE; i++ )
-               {
-                 for( int j = 0; j < 3; j++ )
-                   {
-                     //    fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j]; 
-                     fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i +j -1 ];
-                   }
-               }
-             
-             double diff = 9999;
              int tmpindex = 0;
 
-             for(int k=0; k < 3; k ++)
-               {
-                 //              cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
-                 if(  TMath::Abs(fAmpA[k] - ( maxamp - ped) )  < diff)
-                   {
-                     diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
-                     tmpindex = k; 
-                   }
-               }
-             
+
              Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
              
              if( tmptof < -1 )
@@ -192,124 +169,186 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
              else
                if( tmptof  > -1 && tmptof < 100 )
                  {
-                   tmpindex =1;
+                   tmpindex = 1;
                  }
                else
                  {
                    tmpindex = 2;
                  }
+
              double tof = 0;
-             
-             for(int k=0; k < SAMPLERANGE; k++   )
+           
+             for(int k=0; k < PF::SAMPLERANGE; k++   )
                {
                  tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];   
                }
-             
-             //              cout << __FILE__ << __LINE__ <<  "tofRaw =   "<< tof /  fAmpA[tmpindex]  << endl;
-             
-             // tof = tof /  fAmpA[tmpindex] +  (dt + startbin)*100;
-             
-             if( TMath::Abs(  (maxf - fAmpA[tmpindex])/maxf )  >   0.1 )
+           
+             for( int i=0; i < PF::SAMPLERANGE; i++ )
                {
-                 fAmpA[tmpindex] = maxf;
+                 {
+                   fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i  +tmpindex -1 ];
+                 }
+               }
+             if( TMath::Abs(  (maxf - fAmp  )/maxf )  >   0.1 )
+               {
+                 //      cout << __FILE__ << ":" << __LINE__ << "WARNING: amp was" << fAmp <<", but was changed to "<< maxf << endl;
+                 fAmp = maxf;
                }
-
-             //              timebinOffset
-
-               //      tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex];
-               // tof = ( timebinOffset )*100 - tof/fAmpA[tmpindex]; // ns
-             tof = timebinOffset - 0.01*tof/fAmpA[tmpindex]; // clock ticks
              
-             //      tof = tof/fAmpA[tmpindex];
-
-  
-             return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof, -2, -3 ); 
+             tof = timebinOffset - 0.01*tof/fAmp; // clock ticks
+             // use local-array time for chi2 estimate
+             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 , -5, maxf, -6, -7, -8 ); 
+             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 filenameCoarse[256];
-         char filename[256];
-         
-         int n = j+fNsampleCut;
+         for( int j=0; j < PF::SAMPLERANGE; j++)  
+           {
+             pfv->GetVector( i, j, fPFAmpVectors[i][j] ,  fPFTofVectors[i][j],    
+                             fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j]  ); 
 
-         //      double start = (double)i+0.5;
-         double start = (double)i+0;
+             fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
+                                            fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
+           }
+       }
+    }
+  else
+    {
+      AliFatal( "pfv = ZERO !!!!!!!");
+    } 
+}
 
-         sprintf(filename,        "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
-         sprintf(filenameCoarse,  "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
 
-         FILE *fp  =  fopen(filename, "r");
-         FILE *fpc =  fopen(filenameCoarse, "r");
+void   
+AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
+{
+  //Loading of Peak-Finder  vectors from the 
+  //Offline Condition Database  (OCDB)
+  AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
+  cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
+  entry->PrintMetaData();
+
+  if( entry != 0 )
+    {
+      AliCaloPeakFinderVectors  *pfv = (AliCaloPeakFinderVectors *)entry->GetObject(); 
+      if( pfv == 0 )
+       {
+         cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
+       }
+  
+      CopyVectors( pfv );
+      
+      if( pfv != 0 )
+       {
+         fIsInitialized = true;
+       }
+    }
+}
+
 
-         if( fp == 0 )
+void 
+AliCaloRawAnalyzerPeakFinder::LoadVectorsASCII()
+{
+  //Read in the Peak finder vecors from ASCI files
+  fIsInitialized= true;  
+  const Int_t buffersize = 256;
+  for(int i = 0;  i < PF::MAXSTART ; i++)
+  {
+    for( int j=0; j < PF::SAMPLERANGE; j++)
+    {
+      char filenameCoarse[buffersize];
+      char filename[buffersize];
+      int n = j+fNsampleCut;
+      double start = (double)i+0;
+      
+      snprintf(filename, buffersize,       "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
+      snprintf(filenameCoarse, buffersize, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
+      
+      FILE *fp  =  fopen(filename, "r");
+      FILE *fpc =  fopen(filenameCoarse, "r");
+      
+      if( fp == 0 )
            {
              AliFatal( Form( "could not open file: %s", filename ) );
            }
-       
-         if(fpc == 0)
+      else if(fpc == 0)
            {
              AliFatal( Form( "could not open file: %s", filenameCoarse ) );
            }
-         else
+      else
            {
              for(int m = 0; m < n ; m++ )
-               {
-                 cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
-                 fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
-                 //              fPFAmpVectorsCoarse[i][j][m] = 1;
-                 fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
-               }
-             
+        {
+          fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
+          fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
+        }
              fscanf(fp,   "\n" );
              fscanf(fpc,  "\n" );
-             
              for(int m = 0; m < n ; m++ )
-               {
-                 //  fPFTofVectors[i][j][m] = 1;
-
-                 fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
-                 fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );  
-                 //  fPFTofVectorsCoarse[i][j][m] = 1;  
-               }
+        {
+          fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
+          fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );  
+        }
              
-            
+             fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
+                                      fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
+        
              fclose (fp);
              fclose (fpc);
            }
-       }
     }
+  }
 }
 
 
+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::PolTof( const double rectof ) const
-//
+void 
+AliCaloRawAnalyzerPeakFinder::PrintVectors()
 {
-  static Double_t p0 = -55.69;
-  static Double_t p1 = 3.178;
-  static Double_t p2 = -0.05587;
-  static Double_t p3 = 0.0003185;
-  static Double_t p4 = -7.91E-7;
-  static Double_t p5 = 7.576E-10;
+  for(int i=0; i < 20; i++)
+    {
+      for( int j = 0; j < PF::MAXSTART; j ++ )
+       {
+         for( int k=0; k < PF::SAMPLERANGE; k++ )
+           {
+             cout << fPFAmpVectors[j][k][i] << "\t" ;
+           }
+       }
+      cout << endl;
+    }
+  cout << __FILE__ << ":" << __LINE__ << ":.... DONE !!" << endl;
 }
-*/