EMCAL/AliCaloRawAnalyzerPeakFinder.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by                              *
   3  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
   4  *                                                                        *
   5  * Primary Author: Per Thomas Hille  <perthomas.hille@yale.edu>           *
   6  *                                                                        *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  * Please report bugs to   perthomas.hille@yale.edu                       *
   9  *                                                                        *
  10  * Permission to use, copy, modify and distribute this software and its   *
  11  * documentation strictly for non-commercial purposes is hereby granted   *
  12  * without fee, provided that the above copyright notice appears in all   *
  13  * copies and that both the copyright notice and this permission notice   *
  14  * appear in the supporting documentation. The authors make no claims     *
  15  * about the suitability of this software for any purpose. It is          *
  16  * provided "as is" without express or implied warranty.                  *
  17  **************************************************************************/
  18
  19
  20 // The Peak-Finder algorithm
  21 // The amplitude is extracted  as a
  22 // weighted sum of the samples using the
  23 // best possible weights.
  24 // The wights is calculated only once and the
  25 // Actual extraction of amplitude and peak position
  26 // Is done with a simple vector multiplication, allowing for
  27 // Extreemely fast computations.
  28
  29 #include "AliCaloRawAnalyzerPeakFinder.h"
  30 #include "AliCaloBunchInfo.h"
  31 #include "AliCaloFitResults.h"
  32 #include "TMath.h"
  33 #include "AliLog.h"
  34 #include "AliCDBEntry.h"
  35 #include "AliCDBManager.h"
  36 #include "TFile.h"
  37 #include "AliCaloPeakFinderVectors.h"
  38
  39 using namespace std;
  40
  41
  42 ClassImp( AliCaloRawAnalyzerPeakFinder )
  43
  44
  45 AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),
  46                                                               fAmp(0),
  47                                                               fPeakFinderVectors(0),
  48                                                               fRunOnAlien(false)
  49 {
  50   //Comment
  51   InitOCDB(fRunOnAlien);
  52   fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
  53   ResetVectors();
  54   LoadVectorsOCDB();
  55 }
  56
  57
  58 void
  59 AliCaloRawAnalyzerPeakFinder::InitOCDB(bool alien) const
  60 {
  61   // Setting the default OCDB pathe depending on wether we work locally or on the GRID.
  62   AliCDBManager::Instance()->SetDefaultStorage(  alien == true ? "alien://$ALICE_ROOT/OCDB" : "local://$ALICE_ROOT/OCDB");
  63   AliCDBManager::Instance()->SetRun(100);
  64 }
  65
  66
  67 void
  68 AliCaloRawAnalyzerPeakFinder::ResetVectors()
  69 {
  70   //As name implies
  71   for(int i=0; i < MAXSTART; i++)
  72     {
  73       for(int j=0; j < SAMPLERANGE; j++ )
  74         {
  75           for(int k=0; k < 100; k++ )
  76             {
  77               fPFAmpVectors[i][j][k] = 0;
  78               fPFTofVectors[i][j][k] = 0;
  79               fPFAmpVectorsCoarse[i][j][k] = 0;
  80               fPFTofVectorsCoarse[i][j][k] = 0;
  81             }
  82         }
  83     }
  84 }
  85
  86
  87 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
  88 {
  89   //comment
  90   for(int i=0; i < MAXSTART; i++)
  91     {
  92       for(int j=0; j < SAMPLERANGE; j++ )
  93         {
  94           delete[] fPFAmpVectors[i][j];
  95           delete[] fPFTofVectors[i][j];
  96           delete[] fPFAmpVectorsCoarse[i][j];
  97           delete[] fPFTofVectorsCoarse[i][j];
  98         }
  99     }
 100 }
 101
 102
 103
 104 Double_t
 105 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
 106 {
 107   // Fisrt (coarce) estimate of Amplitude using the Peak-Finder.
 108   // The output of the first iteration is sued to select vectors
 109   // for the second iteration.
 110
 111   Double_t tmpTof = 0;
 112   Double_t tmpAmp= 0;
 113
 114   for(int i=0; i < length; i++)
 115     {
 116       tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
 117       tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
 118     }
 119
 120   tmpTof = tmpTof / tmpAmp ;
 121   return tmpTof;
 122 }
 123
 124
 125 AliCaloFitResults
 126 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
 127 {
 128   // Extracting the amplitude using the Peak-Finder algorithm
 129   // The amplitude is a weighted sum of the samples using
 130   // optimum weights.
 131
 132   short maxampindex; //index of maximum amplitude
 133   short maxamp; //Maximum amplitude
 134   fAmp = 0;
 135   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
 136
 137   if( index >= 0)
 138     {
 139       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
 140       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
 141       short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1);
 142
 143       if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
 144         {
 145           return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
 146         }
 147       else if ( maxf >= fAmpCut )
 148         {
 149           int first = 0;
 150           int last = 0;
 151           short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
 152           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last);
 153           int nsamples =  last - first;
 154
 155           if( ( nsamples  )  >= fNsampleCut ) // no if statement needed really; keep for readability
 156             {
 157               int startbin = bunchvector.at(index).GetStartBin();
 158               int n = last - first;
 159               int pfindex = n - fNsampleCut;
 160               pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
 161
 162               int dt =  maxampindex - startbin -2;
 163               int tmpindex = 0;
 164
 165               Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
 166
 167               if( tmptof < -1 )
 168                 {
 169                   tmpindex = 0;
 170                 }
 171               else
 172                 if( tmptof  > -1 && tmptof < 100 )
 173                   {
 174                     tmpindex = 1;
 175                   }
 176                 else
 177                   {
 178                     tmpindex = 2;
 179                   }
 180
 181               double tof = 0;
 182
 183               for(int k=0; k < SAMPLERANGE; k++   )
 184                 {
 185                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];
 186                 }
 187
 188               for( int i=0; i < SAMPLERANGE; i++ )
 189                 {
 190                   {
 191                     fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i  +tmpindex -1 ];
 192                   }
 193                 }
 194               if( TMath::Abs(  (maxf - fAmp  )/maxf )  >   0.1 )
 195                 {
 196                   fAmp = maxf;
 197                 }
 198
 199               tof = timebinOffset - 0.01*tof/fAmp; // clock ticks
 200
 201               // use local-array time for chi2 estimate
 202               Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
 203               Int_t ndf = last - first - 1; // nsamples - 2
 204               return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fAmp, tof,
 205                                         timebinOffset, chi2, ndf,
 206                                         AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 207             }
 208           else
 209             {
 210               Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
 211               Int_t ndf = last - first - 1; // nsamples - 2
 212               return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kCrude, maxf, timebinOffset,
 213                                         timebinOffset, chi2, ndf, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 214             }
 215         } // ampcut
 216     }
 217   return  AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
 218 }
 219
 220
 221 void
 222 AliCaloRawAnalyzerPeakFinder::CopyVectors(const AliCaloPeakFinderVectors *const pfv )
 223 {
 224   // As name implies
 225
 226   if ( pfv != 0)
 227     {
 228       for(int i = 0;  i < MAXSTART ; i++)
 229         {
 230           for( int j=0; j < SAMPLERANGE; j++)
 231             {
 232               pfv->GetVector( i, j, fPFAmpVectors[i][j] ,  fPFTofVectors[i][j],
 233                               fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j]  );
 234               fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
 235                                              fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
 236             }
 237         }
 238     }
 239   else
 240     {
 241       AliFatal( "pfv = ZERO !!!!!!!");
 242     }
 243 }
 244
 245
 246
 247 void
 248 AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
 249 {
 250   //Loading of Peak-Finder  vectors from the
 251   //Offline Condition Database  (OCDB)
 252   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
 253
 254   if( entry != 0 )
 255     {
 256       AliCaloPeakFinderVectors  *pfv = (AliCaloPeakFinderVectors *)entry->GetObject();
 257       CopyVectors( pfv );
 258      }
 259 }
 260
 261
 262 void
 263 AliCaloRawAnalyzerPeakFinder::LoadVectorsASCII()
 264 {
 265   //Read in the Peak finder vecors from ASCI files
 266   for(int i = 0;  i < MAXSTART ; i++)
 267     {
 268       for( int j=0; j < SAMPLERANGE; j++)
 269         {
 270           char filenameCoarse[256];
 271           char filename[256];
 272           int n = j+fNsampleCut;
 273           double start = (double)i+0;
 274
 275           sprintf(filename,        "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
 276           sprintf(filenameCoarse,  "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
 277
 278           FILE *fp  =  fopen(filename, "r");
 279           FILE *fpc =  fopen(filenameCoarse, "r");
 280
 281           if( fp == 0 )
 282             {
 283               AliFatal( Form( "could not open file: %s", filename ) );
 284             }
 285
 286           if(fpc == 0)
 287             {
 288               AliFatal( Form( "could not open file: %s", filenameCoarse ) );
 289             }
 290           else
 291             {
 292               for(int m = 0; m < n ; m++ )
 293                 {
 294                   fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
 295                   fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
 296                 }
 297               fscanf(fp,   "\n" );
 298               fscanf(fpc,  "\n" );
 299               for(int m = 0; m < n ; m++ )
 300                 {
 301                   fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
 302                   fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );
 303                 }
 304
 305               fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
 306                                              fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
 307
 308               fclose (fp);
 309               fclose (fpc);
 310             }
 311         }
 312     }
 313 }
 314
 315
 316 void
 317 AliCaloRawAnalyzerPeakFinder::WriteRootFile() const
 318 {
 319   // Utility function to write Peak-Finder vectors to an root file
 320   // The output is used to create an OCDB entry.
 321   fPeakFinderVectors->PrintVectors();
 322   TFile *f = new TFile("peakfindervectors2.root",  "recreate" );
 323   fPeakFinderVectors->Write();
 324   f->Close();
 325   delete f;
 326 }