EMCAL/AliCaloRawAnalyzerPeakFinder.cxx

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