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