EMCAL/EMCALraw/AliCaloRawAnalyzer.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 p.t.hille@fys.uio.no                             *
  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
  21 // Base class for extraction
  22 // of signal amplitude and peak position
  23 // From CALO Calorimeter RAW data (from the RCU)
  24 // Contains some utilities for preparing / selecting
  25 // Signals suitable for signal extraction
  26 // By derived classes
  27
  28 #include "AliLog.h"
  29 #include "AliCaloRawAnalyzer.h"
  30 #include "AliCaloBunchInfo.h"
  31 #include "AliCaloFitResults.h"
  32 #include "TMath.h"
  33 #include <iostream>
  34 using namespace std;
  35
  36 ClassImp(AliCaloRawAnalyzer)
  37
  38 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) :  TObject(),
  39   fMinTimeIndex(-1),
  40   fMaxTimeIndex(-1),
  41   fFitArrayCut(5),
  42   fAmpCut(4),
  43   fNsampleCut(5),
  44   fOverflowCut(950),
  45   fNsamplePed(3),
  46   fIsZerosupressed( false ),
  47   fVerbose( false ),
  48   fAlgo(Algo::kNONE),
  49   fL1Phase(0),
  50   fAmp(0),
  51   fTof(0),
  52   fTau( EMCAL::TAU ),
  53   fFixTau( true )
  54 {
  55   // Ctor
  56
  57   snprintf(fName,     256, "%s", name);
  58   snprintf(fNameShort,256, "%s", nameshort);
  59
  60   for(int i=0; i < ALTROMAXSAMPLES; i++ )
  61     {
  62       fReversed[i] = 0;
  63     }
  64 }
  65
  66
  67 void
  68 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
  69 {
  70   //Require that the bin if the maximum ADC value is between min and max (timebin)
  71
  72   if(  ( min > max ) || min > ALTROMAXSAMPLES  || max > ALTROMAXSAMPLES  )
  73     {
  74       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) );
  75     }
  76   else
  77     {
  78       fMinTimeIndex  = min;
  79       fMaxTimeIndex  = max;
  80     }
  81 }
  82
  83
  84 UShort_t
  85 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
  86 {
  87   // Get maximum of array
  88
  89   UShort_t tmpmax  = data[0];
  90
  91   for(int i=0; i < length; i++)
  92     {
  93       if( tmpmax  <  data[i] )
  94         {
  95           tmpmax = data[i];
  96         }
  97     }
  98   return tmpmax;
  99 }
 100
 101
 102 void
 103 AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, int length, short maxindex,
 104                                     int * first, int * last, int cut) const
 105 {
 106   //Selection of subset of data from one bunch that will be used for fitting or
 107   //Peak finding.  Go to the left and right of index of the maximum time bin
 108   //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
 109   int tmpfirst =  maxindex;
 110   int tmplast  =  maxindex;
 111   Double_t prevFirst =  data[maxindex];
 112   Double_t prevLast  =  data[maxindex];
 113   bool firstJump = false;
 114   bool lastJump = false;
 115
 116   while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) )
 117     {
 118       // jump check:
 119       if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
 120         if ( data[tmpfirst] >= prevFirst) {
 121           firstJump = true;
 122         }
 123       }
 124       prevFirst = data[tmpfirst];
 125       tmpfirst -- ;
 126     }
 127
 128   while( (tmplast < length) && (data[tmplast] >=  cut ) && (!lastJump) )
 129     {
 130       // jump check:
 131       if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
 132         if ( data[tmplast] >= prevLast) {
 133           lastJump = true;
 134         }
 135       }
 136       prevLast = data[tmplast];
 137       tmplast ++;
 138     }
 139
 140   // we keep one pre- or post- sample if we can (as in online)
 141   // check though if we ended up on a 'jump', or out of bounds: if so, back up
 142   if (firstJump || tmpfirst<0) tmpfirst ++;
 143   if (lastJump || tmplast>=length) tmplast --;
 144
 145   *first = tmpfirst;
 146   *last =  tmplast;
 147
 148   return;
 149 }
 150
 151
 152
 153 Float_t
 154 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch,
 155                                            UInt_t /*altrocfg1*/, UInt_t /*altrocfg2*/,
 156                                            double *outarray ) const
 157 {
 158   //Time sample comes in reversed order, revers them back
 159   //Subtract the baseline based on content of altrocfg1 and altrocfg2.
 160
 161   Int_t length =  bunch->GetLength();
 162   const UShort_t *sig = bunch->GetData();
 163
 164   double ped = EvaluatePedestal( sig , length);
 165
 166   for( int i=0; i < length; i++ )
 167     {
 168       outarray[i] = sig[length -i -1] - ped;
 169     }
 170
 171   return ped;
 172 }
 173
 174
 175
 176 Float_t
 177 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * data, int /*length*/ ) const
 178 {
 179   //  Pedestal evaluation if not zero suppressed
 180
 181   double tmp = 0;
 182
 183   if( fIsZerosupressed == false )
 184     {
 185       for(int i=0; i < fNsamplePed; i++ )
 186         {
 187           tmp += data[i];
 188         }
 189    }
 190
 191   return  tmp/fNsamplePed;
 192  }
 193
 194
 195 short
 196 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo * bunch , int * maxindex ) const
 197 {
 198   // Get maximum in bunch array
 199
 200   short tmpmax = -1;
 201   int tmpindex = -1;
 202   const UShort_t *sig = bunch->GetData();
 203
 204   for(int i=0; i < bunch->GetLength(); i++ )
 205     {
 206       if( sig[i] > tmpmax )
 207         {
 208           tmpmax   =  sig[i];
 209           tmpindex =  i;
 210         }
 211     }
 212
 213   if(maxindex != 0 )
 214     {
 215       //   *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
 216        *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
 217     }
 218
 219   return  tmpmax;
 220 }
 221
 222
 223 bool
 224 AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
 225 {
 226   // a bunch is considered invalid if the maximum is in the first or last time-bin
 227   short tmpmax = -1;
 228   int tmpindex = -1;
 229   const UShort_t *sig = bunch->GetData();
 230
 231   for(int i=0; i < bunch->GetLength(); i++ )
 232     {
 233       if( sig[i] > tmpmax )
 234         {
 235           tmpmax   =  sig[i];
 236           tmpindex =  i;
 237         }
 238     }
 239
 240   bool bunchOK = true;
 241   if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
 242     {
 243       bunchOK = false;
 244     }
 245
 246   return  bunchOK;
 247 }
 248
 249
 250 int
 251 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,
 252                                  short * maxampbin, short * maxamplitude )
 253 {
 254   //We select the bunch with the highest amplitude unless any time constraints is set
 255   short max = -1;
 256   short bunchindex = -1;
 257   short maxall = -1;
 258   int indx = -1;
 259
 260   for(unsigned int i=0; i < bunchvector.size(); i++ )
 261     {
 262       max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
 263       if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
 264         {
 265           if( max > maxall )
 266             {
 267               maxall = max;
 268               bunchindex = i;
 269               *maxampbin     = indx;
 270               *maxamplitude  = max;
 271             }
 272         }
 273     }
 274
 275   if (bunchindex >= 0) {
 276     bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
 277     if (! bunchOK) {
 278       bunchindex = -1;
 279     }
 280   }
 281
 282   return  bunchindex;
 283 }
 284
 285
 286 bool
 287 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
 288 {
 289   // Ckeck if the index of the max ADC vaue is consistent with trigger.
 290   if( ( mintindx  < 0 && maxtindx   < 0) ||maxtindx  < 0 )
 291     {
 292       return true;
 293     }
 294   return ( maxindex < maxtindx ) && ( maxindex > mintindx  ) ? true : false;
 295 }
 296
 297
 298 void
 299 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr )
 300 {
 301   // Print bunch vector info
 302   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
 303   cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
 304
 305   for(unsigned int i=0; i < bvctr.size() ; i++ )
 306     {
 307       PrintBunch( bvctr.at(i) );
 308       cout << " bunch = "  <<  i  << endl;
 309     }
 310   cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
 311 }
 312
 313
 314 void
 315 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch )
 316 {
 317   // Print bunch infor
 318   cout << __FILE__ << ":" << __LINE__ << endl;
 319   cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
 320   const UShort_t *sig =  bunch.GetData();
 321
 322   for ( Int_t j = 0; j <  bunch.GetLength();  j++)
 323     {
 324       printf("%d\t", sig[j] );
 325     }
 326   cout << endl;
 327 }
 328
 329
 330 Double_t
 331 AliCaloRawAnalyzer::CalculateChi2( Double_t amp,    Double_t time,
 332                                    Int_t    first,  Int_t    last,
 333                                    Double_t adcErr, Double_t tau) const
 334 {
 335   //   Input:
 336   //   amp   - max amplitude;
 337   //   time    - time of max amplitude;
 338   //   first, last - sample array indices to be used
 339   //   adcErr   - nominal error of amplitude measurement (one value for all channels)
 340   //           if adcErr<0 that mean adcErr=1.
 341   //   tau   - filter time response (in timebin units)
 342   // Output:
 343   //   chi2 - chi2
 344
 345   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
 346     // or, first is negative, the indices are not valid
 347     return Ret::kDummy;
 348   }
 349
 350   int nsamples =  last - first + 1;
 351   // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time);
 352
 353   Int_t x = 0;
 354   Double_t chi2 = 0;
 355   Double_t dy = 0.0, xx = 0.0, f=0.0;
 356
 357   for (Int_t i=0; i<nsamples; i++) {
 358     x     = first + i; // timebin
 359     xx    = (x - time + tau) / tau; // help variable
 360     f     = 0;
 361     if (xx > 0) {
 362       f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
 363     }
 364     dy    = fReversed[x] - f;
 365     chi2 += dy*dy;
 366     // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
 367   }
 368
 369   if (adcErr>0.0) { // weight chi2
 370     chi2 /= (adcErr*adcErr);
 371   }
 372   return chi2;
 373 }
 374
 375
 376 void
 377 AliCaloRawAnalyzer::CalculateMeanAndRMS(Int_t first, Int_t last,
 378                                         Double_t & mean, Double_t & rms)
 379 {
 380   //   Input:
 381   //   first, last - sample array indices to be used
 382   // Output:
 383   //   mean and RMS of samples
 384   //
 385   // To possibly be used to differentiate good signals from bad before fitting
 386   //
 387   mean = Ret::kDummy;
 388   rms =  Ret::kDummy;
 389
 390   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
 391     // or, first is negative, the indices are not valid
 392     return;
 393   }
 394
 395   int nsamples =  last - first + 1;
 396   //  printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
 397
 398   int x = 0;
 399   Double_t sampleSum = 0; // sum of samples
 400   Double_t squaredSampleSum = 0; // sum of samples squared
 401
 402   for (Int_t i=0; i<nsamples; i++) {
 403     x = first + i;
 404     sampleSum += fReversed[x];
 405     squaredSampleSum += (fReversed[x] * fReversed[x]);
 406   }
 407
 408   mean = sampleSum / nsamples;
 409   Double_t squaredMean = squaredSampleSum / nsamples;
 410   // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
 411   rms = sqrt(squaredMean - mean*mean);
 412
 413   return;
 414 }
 415
 416
 417
 418 int
 419 AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector,
 420                                           UInt_t altrocfg1, UInt_t altrocfg2, Int_t & index,
 421                                           Float_t & maxf, short & maxamp,
 422                                           short & maxrev, Float_t & ped,
 423                                           int & first, int & last, int acut )
 424 {
 425   // method to do the selection of what should possibly be fitted
 426
 427   int nsamples  = 0;
 428   short maxampindex = 0;
 429   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
 430
 431
 432   if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
 433     {
 434       // use more convenient numbering and possibly subtract pedestal
 435       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
 436       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
 437
 438       if ( maxf >= acut  ) // possibly significant signal
 439         {
 440           // select array around max to possibly be used in fit
 441           maxrev = maxampindex - bunchvector.at(index).GetStartBin();
 442           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, acut );
 443
 444           // sanity check: maximum should not be in first or last bin
 445           // if we should do a fit
 446           if (first!=maxrev && last!=maxrev)
 447             {
 448               // calculate how many samples we have
 449               nsamples =  last - first + 1;
 450             }
 451         }
 452     }
 453
 454   return nsamples;
 455 }
 456