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                                                                                    fNsamplePed(3),
  44                                                                                    fIsZerosupressed( false ),
  45                                                                                    fVerbose( false )
  46 {
  47   //Comment
  48   sprintf(fName, "%s", name);
  49   sprintf(fNameShort, "%s", nameshort);
  50
  51   for(int i=0; i < MAXSAMPLES; i++ )
  52     {
  53       fReversed[i] = 0;
  54     }
  55 }
  56
  57 AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
  58 {
  59
  60 }
  61
  62
  63 void
  64 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
  65 {
  66   //Require that the bin if the maximum ADC value is between min and max (timebin)
  67   if(  ( min > max ) || min > MAXSAMPLES  || max > MAXSAMPLES  )
  68     {
  69       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) );
  70     }
  71   else
  72     {
  73       fMinTimeIndex  = min;
  74       fMaxTimeIndex  = max;
  75     }
  76 }
  77
  78
  79 UShort_t
  80 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
  81 {
  82   //------------
  83   UShort_t tmpmax  = data[0];
  84
  85   for(int i=0; i < length; i++)
  86     {
  87       if( tmpmax  <  data[i] )
  88         {
  89           tmpmax = data[i];
  90         }
  91     }
  92   return tmpmax;
  93 }
  94
  95
  96 void
  97 AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first,  int *const last ) const
  98 {
  99   //Selection of subset of data from one bunch that will be used for fitting or
 100   //Peak finding.  Go to the left and right of index of the maximum time bin
 101   //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
 102   int tmpfirst =  maxindex;
 103   int tmplast  =  maxindex;
 104   Double_t valFirst =  fData[maxindex];
 105   Double_t valLast  =  fData[maxindex];
 106
 107   while( (tmpfirst >= 0) && (fData[tmpfirst] >= fFitArrayCut) &&
 108          (fData[tmpfirst]<valFirst || tmpfirst==maxindex) )
 109     {
 110       valFirst = fData[tmpfirst];
 111       tmpfirst -- ;
 112     }
 113
 114   while( (tmplast < length) && (fData[tmplast] >= fFitArrayCut) &&
 115          (fData[tmplast]<valLast || tmplast==maxindex) )
 116     {
 117       valLast = fData[tmplast];
 118       tmplast ++;
 119     }
 120
 121   *first = tmpfirst +1;
 122   *last =  tmplast -1;
 123   return;
 124 }
 125
 126
 127
 128 Float_t
 129 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/, double *outarray ) const
 130 {
 131   //Time sample comes in reversed order, revers them back
 132   //Subtract the baseline based on content of altrocfg1 and altrocfg2.
 133   Int_t length =  bunch->GetLength();
 134   const UShort_t *sig = bunch->GetData();
 135
 136   double ped = EvaluatePedestal( sig , length);
 137
 138   for( int i=0; i < length; i++ )
 139     {
 140       outarray[i] = sig[length -i -1] - ped;
 141     }
 142
 143   return ped;
 144 }
 145
 146
 147
 148 Float_t
 149 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
 150 {
 151   //  double ped = 0;
 152   double tmp = 0;
 153
 154   if( fIsZerosupressed == false )
 155     {
 156       for(int i=0; i < fNsamplePed; i++ )
 157         {
 158           tmp += data[i];
 159         }
 160    }
 161
 162   return  tmp/fNsamplePed;
 163  }
 164
 165
 166 short
 167 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const
 168 {
 169   //comment
 170   short tmpmax = -1;
 171   int tmpindex = -1;
 172   const UShort_t *sig = bunch->GetData();
 173
 174   for(int i=0; i < bunch->GetLength(); i++ )
 175     {
 176       if( sig[i] > tmpmax )
 177         {
 178           tmpmax   =  sig[i];
 179           tmpindex =  i;
 180         }
 181     }
 182
 183   if(maxindex != 0 )
 184     {
 185       //   *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
 186        *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
 187     }
 188
 189   return  tmpmax;
 190 }
 191
 192
 193 bool
 194 AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
 195 {
 196   // a bunch is considered invalid if the maximum is in the first or last time-bin
 197   short tmpmax = -1;
 198   int tmpindex = -1;
 199   const UShort_t *sig = bunch->GetData();
 200
 201   for(int i=0; i < bunch->GetLength(); i++ )
 202     {
 203       if( sig[i] > tmpmax )
 204         {
 205           tmpmax   =  sig[i];
 206           tmpindex =  i;
 207         }
 208     }
 209
 210   bool bunchOK = true;
 211   if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
 212     {
 213       bunchOK = false;
 214     }
 215
 216   return  bunchOK;
 217 }
 218
 219
 220 int
 221 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
 222 {
 223   //We select the bunch with the highest amplitude unless any time constraints is set
 224   short max = -1;
 225   short bunchindex = -1;
 226   short maxall = -1;
 227   int indx = -1;
 228
 229   for(unsigned int i=0; i < bunchvector.size(); i++ )
 230     {
 231       max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
 232       if( IsInTimeRange( indx) )
 233         {
 234           if( max > maxall )
 235             {
 236               maxall = max;
 237               bunchindex = i;
 238               *maxampbin     = indx;
 239               *maxamplitude  = max;
 240             }
 241         }
 242     }
 243
 244   if (bunchindex >= 0) {
 245     bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
 246     if (! bunchOK) {
 247       bunchindex = -1;
 248     }
 249   }
 250
 251   return  bunchindex;
 252 }
 253
 254
 255 bool
 256 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex  ) const
 257 {
 258   // Ckeck if the index of the max ADC vaue is consistent with trigger.
 259   if( ( fMinTimeIndex  < 0 && fMaxTimeIndex  < 0) ||fMaxTimeIndex  < 0 )
 260     {
 261       return true;
 262     }
 263   return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex  ) ? true : false;
 264 }
 265
 266
 267 void
 268 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
 269 {
 270   //comment
 271   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
 272   cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
 273
 274   for(unsigned int i=0; i < bvctr.size() ; i++ )
 275     {
 276       PrintBunch( bvctr.at(i) );
 277       cout << " bunch = "  <<  i  << endl;
 278     }
 279   cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
 280 }
 281
 282
 283 void
 284 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
 285 {
 286   //comment
 287   cout << __FILE__ << ":" << __LINE__ << endl;
 288   cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
 289   const UShort_t *sig =  bunch.GetData();
 290
 291   for ( Int_t j = 0; j <  bunch.GetLength();  j++)
 292     {
 293       printf("%d\t", sig[j] );
 294     }
 295   cout << endl;
 296 }
 297
 298
 299 AliCaloFitResults
 300 AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo>  &/*bunchvector*/, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
 301 { // method to do the selection of what should possibly be fitted
 302   // not implemented for base class
 303   return AliCaloFitResults( 0, 0 );
 304 }
 305
 306
 307 int
 308 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)
 309 { // method to do the selection of what should possibly be fitted
 310   int nsamples  = 0;
 311   short maxampindex = 0;
 312   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
 313
 314
 315   if( index >= 0 && maxamp > fAmpCut) // something valid was found, and non-zero amplitude
 316     {
 317       // use more convenient numbering and possibly subtract pedestal
 318       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
 319       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
 320
 321       if ( maxf > fAmpCut ) // possibly significant signal
 322         {
 323           // select array around max to possibly be used in fit
 324           maxrev = maxampindex - bunchvector.at(index).GetStartBin();
 325           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
 326
 327           // sanity check: maximum should not be in first or last bin
 328           // if we should do a fit
 329           if (first!=maxrev && last!=maxrev) {
 330             // calculate how many samples we have
 331             nsamples =  last - first + 1;
 332           }
 333         }
 334     }
 335
 336   return nsamples;
 337 }
 338