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