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