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