EMCAL/AliCaloRawAnalyzerKStandard.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 <p.t.hille@fys.uio.no>                *
   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 // Extraction of amplitude and peak position
  21 // FRom CALO raw data using
  22 // least square fit for the
  23 // Moment assuming identical and
  24 // independent errors (equivalent with chi square)
  25 //
  26
  27 #include "AliCaloRawAnalyzerKStandard.h"
  28 #include "AliCaloBunchInfo.h"
  29 #include "AliCaloFitResults.h"
  30 #include "AliLog.h"
  31 #include "TMath.h"
  32 #include <stdexcept>
  33 #include <iostream>
  34 #include "TF1.h"
  35 #include "TGraph.h"
  36 #include "TRandom.h"
  37
  38
  39 using namespace std;
  40
  41
  42 #define BAD 4  //CRAP PTH
  43
  44 ClassImp( AliCaloRawAnalyzerKStandard )
  45
  46
  47 AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzer("Chi Square ( kStandard )", "KStandard"),
  48                                                  fkEulerSquared(7.389056098930650227),
  49                                                  fTf1(0),
  50                                                  fTau(2.35),
  51                                                  fFixTau(kTRUE)
  52 {
  53
  54   fAlgo = Algo::kStandard;
  55   //comment
  56   for(int i=0; i < ALTROMAXSAMPLES; i++)
  57     {
  58       fXaxis[i] = i;
  59     }
  60
  61   fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])",  0, 30 );
  62   if (fFixTau)
  63     {
  64     fTf1->FixParameter(2, fTau);
  65     }
  66   else
  67     {
  68       fTf1->ReleaseParameter(2); // allow par. to vary
  69       fTf1->SetParameter(2, fTau);
  70     }
  71 }
  72
  73
  74 AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
  75 {
  76   delete fTf1;
  77 }
  78
  79
  80
  81 AliCaloFitResults
  82 AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo>  &bunchlist, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  83 {
  84
  85   Float_t pedEstimate  = 0;
  86   short maxADC = 0;
  87   Int_t first = 0;
  88   Int_t last = 0;
  89   Int_t bunchIndex = 0;
  90   Float_t ampEstimate = 0;
  91   short timeEstimate  = 0;
  92   Float_t time = 0;
  93   Float_t amp=0;
  94   Float_t chi2 = 0;
  95   Int_t ndf = 0;
  96   Bool_t fitDone = kFALSE;
  97
  98
  99   int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
 100                                         maxADC, timeEstimate, pedEstimate, first, last,   fAmpCut );
 101
 102
 103   if (ampEstimate >= fAmpCut  )
 104     {
 105       time = timeEstimate;
 106       Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
 107       amp = ampEstimate;
 108
 109       if ( nsamples > 1 && maxADC< OVERFLOWCUT )
 110         {
 111           FitRaw(first, last, amp, time, chi2, fitDone);
 112           time += timebinOffset;
 113           timeEstimate += timebinOffset;
 114           ndf = nsamples - 2;
 115         }
 116     }
 117   if ( fitDone )
 118     {
 119       Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
 120       Float_t timeDiff = time - timeEstimate;
 121
 122       if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
 123         {
 124           amp     = ampEstimate;
 125           time    = timeEstimate;
 126           fitDone = kFALSE;
 127         }
 128     }
 129   if (amp >= fAmpCut )
 130     {
 131       if ( ! fitDone)
 132         {
 133           amp += (0.5 - gRandom->Rndm());
 134         }
 135       //Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 136       // lowGain  = in.IsLowGain();
 137
 138       time = time * TIMEBINWITH;
 139
 140       /////////////!!!!!!!!!time -= in.GetL1Phase();
 141
 142       time -= fL1Phase;
 143
 144       //    AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
 145       //  AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf);
 146
 147
 148       return AliCaloFitResults( -99, -99, fAlgo , amp, time,
 149                                 time, chi2, ndf, Ret::kDummy );
 150
 151
 152       //        AliCaloFitSubarray(index, maxrev, first, last));
 153
 154     }
 155
 156
 157   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 158 }
 159
 160
 161
 162
 163 /*
 164    return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
 165                                         timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 166             }
 167         } // ampcut
 168     }
 169   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 170 */
 171
 172
 173 /*
 174   // Extracting signal parameters using fitting
 175   short maxampindex; //index of maximum amplitude
 176   short maxamp; //Maximum amplitude
 177   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
 178
 179   if( index >= 0)
 180     {
 181       Float_t  ped  = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
 182       Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
 183       short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
 184       // timebinOffset is timebin value at maximum (maxrev)
 185       short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
 186       if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
 187         {
 188           return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
 189         }
 190       else if ( maxf >= fAmpCut )
 191         {
 192           int first = 0;
 193           int last = 0;
 194           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, fFitArrayCut);
 195           int nsamples =  last - first + 1;
 196
 197           if( ( nsamples  )  >= fNsampleCut )
 198             {
 199               Float_t tmax = (maxrev - first); // local tmax estimate
 200               TGraph *graph =  new TGraph(  nsamples, fXaxis,  &fReversed[first] );
 201               fTf1->SetParameter(0, maxf*fkEulerSquared );
 202               fTf1->SetParameter(1, tmax - fTau);
 203               // set rather loose parameter limits
 204               fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
 205               fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
 206
 207               if (fFixTau) {
 208                 fTf1->FixParameter(2, fTau);
 209               }
 210               else {
 211                 fTf1->ReleaseParameter(2); // allow par. to vary
 212                 fTf1->SetParameter(2, fTau);
 213               }
 214
 215               Short_t tmpStatus = 0;
 216               try {
 217                 tmpStatus =  graph->Fit(fTf1, "Q0RW");
 218               }
 219               catch (const std::exception & e) {
 220                 AliError( Form("TGraph Fit exception %s", e.what()) );
 221                 return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
 222                                           timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 223               }
 224
 225               if( fVerbose == true )
 226                 {
 227                   AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
 228                   PrintFitResult( fTf1 ) ;
 229                 }
 230               // global tmax
 231               tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
 232                 + fTf1->GetParameter(2); // +tau, makes sum tmax
 233
 234                 delete graph;
 235                 return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
 236                                           fTf1->GetParameter(0)/fkEulerSquared,
 237                                           tmax,
 238                                           timebinOffset,
 239                                           fTf1->GetChisquare(),
 240                                           fTf1->GetNDF(),
 241                                           Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 242
 243                 //     delete graph;
 244
 245             }
 246           else
 247             {
 248
 249               Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
 250               Int_t ndf = last - first - 1; // nsamples - 2
 251               return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
 252                                         timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 253             }
 254         } // ampcut
 255     }
 256   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 257
 258 }
 259 */
 260
 261
 262
 263 void
 264 AliCaloRawAnalyzerKStandard::PrintFitResult(const TF1 *f) const
 265 {
 266   //comment
 267   cout << endl;
 268   cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
 269   cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
 270   cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
 271   cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
 272   //  cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus()  << ",.. !!!!" << endl << endl;
 273   cout << endl << endl;
 274 }
 275
 276
 277
 278
 279
 280 //____________________________________________________________________________
 281 void
 282  AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
 283 { // Fits the raw signal time distribution
 284
 285   //--------------------------------------------------
 286   //Do the fit, different fitting algorithms available
 287   //--------------------------------------------------
 288
 289   // fprintf(fp, "%s:%d:%s\n", __FILE__, __LINE__, __FUNCTION__ );
 290
 291   int nsamples = lastTimeBin - firstTimeBin + 1;
 292   fitDone = kFALSE;
 293
 294   // switch(fFittingAlgorithm)
 295   //   {
 296   //  case Algo::kStandard:
 297   //    {
 298   if (nsamples < 3) { return; } // nothing much to fit
 299   //printf("Standard fitter \n");
 300
 301   // Create Graph to hold data we will fit
 302
 303   TGraph *gSig =  new TGraph( nsamples);
 304
 305   for (int i=0; i<nsamples; i++)
 306     {
 307       Int_t timebin = firstTimeBin + i;
 308       gSig->SetPoint(i, timebin, GetReversed(timebin));
 309     }
 310
 311   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
 312   signalF->SetParameters(10.,5., TAU  ,ORDER,0.); //set all defaults once, just to be safe
 313   signalF->SetParNames("amp","t0","tau","N","ped");
 314   signalF->FixParameter(2,TAU); // tau in units of time bin
 315   signalF->FixParameter(3,ORDER); // order
 316   signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here
 317   signalF->SetParameter(1, time);
 318   signalF->SetParameter(0, amp);
 319   // set rather loose parameter limits
 320   signalF->SetParLimits(0, 0.5*amp, 2*amp );
 321   signalF->SetParLimits(1, time - 4, time + 4);
 322
 323   try {
 324     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
 325     // assign fit results
 326     amp  = signalF->GetParameter(0);
 327     time = signalF->GetParameter(1);
 328     chi2 = signalF->GetChisquare();
 329     fitDone = kTRUE;
 330   }
 331   catch (const std::exception & e) {
 332     AliError( Form("TGraph Fit exception %s", e.what()) );
 333     // stay with default amp and time in case of exception, i.e. no special action required
 334     fitDone = kFALSE;
 335   }
 336   delete signalF;
 337
 338   //printf("Std   : Amp %f, time %g\n",amp, time);
 339   delete gSig; // delete TGraph
 340
 341   //    break;
 342   //    }//kStandard Fitter
 343   //----------------------------
 344
 345   /*
 346     case Algo::kLogFit:
 347     {
 348     if (nsamples < 3) { return; } // nothing much to fit
 349     //printf("LogFit \n");
 350
 351     // Create Graph to hold data we will fit
 352     TGraph *gSigLog =  new TGraph( nsamples);
 353     for (int i=0; i<nsamples; i++) {
 354     Int_t timebin = firstTimeBin + i;
 355     gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) );
 356     }
 357
 358     TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0,  TIMEBINS , 5);
 359     signalFLog->SetParameters(2.3, 5.,TAU,ORDER,0.); //set all defaults once, just to be safe
 360     signalFLog->SetParNames("amplog","t0","tau","N","ped");
 361     signalFLog->FixParameter(2,TAU); // tau in units of time bin
 362     signalFLog->FixParameter(3, ORDER); // order
 363     signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here
 364     signalFLog->SetParameter(1, time);
 365     if (amp>=1) {
 366     signalFLog->SetParameter(0, TMath::Log(amp));
 367     }
 368
 369     gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
 370
 371     // assign fit results
 372     Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp
 373     amp = TMath::Exp(amplog);
 374     time = signalFLog->GetParameter(1);
 375     fitDone = kTRUE;
 376
 377     delete signalFLog;
 378     //printf("LogFit: Amp %f, time %g\n",amp, time);
 379     delete gSigLog;
 380     break;
 381     } //kLogFit
 382       //----------------------------
 383       //----------------------------
 384       }//switch fitting algorithms
 385   */
 386   return;
 387 }
 388
 389
 390 //__________________________________________________________________
 391 void
 392 AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const
 393 {
 394   //BEG YS alternative methods to calculate the amplitude
 395   Double_t * ymx = gSig->GetX() ;
 396   Double_t * ymy = gSig->GetY() ;
 397   const Int_t kN = 3 ;
 398   Double_t ymMaxX[kN] = {0., 0., 0.} ;
 399   Double_t ymMaxY[kN] = {0., 0., 0.} ;
 400   Double_t ymax = 0. ;
 401   // find the maximum amplitude
 402   Int_t ymiMax = 0 ;
 403   for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
 404     if (ymy[ymi] > ymMaxY[0] ) {
 405       ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
 406       ymMaxX[0] = ymx[ymi] ;
 407       ymiMax = ymi ;
 408     }
 409   }
 410   // find the maximum by fitting a parabola through the max and the two adjacent samples
 411   if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
 412     ymMaxY[1] = ymy[ymiMax+1] ;
 413     ymMaxY[2] = ymy[ymiMax-1] ;
 414     ymMaxX[1] = ymx[ymiMax+1] ;
 415     ymMaxX[2] = ymx[ymiMax-1] ;
 416     if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
 417       //fit a parabola through the 3 points y= a+bx+x*x*x
 418       Double_t sy = 0 ;
 419       Double_t sx = 0 ;
 420       Double_t sx2 = 0 ;
 421       Double_t sx3 = 0 ;
 422       Double_t sx4 = 0 ;
 423       Double_t sxy = 0 ;
 424       Double_t sx2y = 0 ;
 425       for (Int_t i = 0; i < kN ; i++) {
 426         sy += ymMaxY[i] ;
 427         sx += ymMaxX[i] ;
 428         sx2 += ymMaxX[i]*ymMaxX[i] ;
 429         sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 430         sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 431         sxy += ymMaxX[i]*ymMaxY[i] ;
 432         sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
 433       }
 434       Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
 435       Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
 436       Double_t c  = cN / cD ;
 437       Double_t b  = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
 438       Double_t a  = (sy-b*sx-c*sx2)/kN  ;
 439       Double_t xmax = -b/(2*c) ;
 440       ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
 441       amp = ymax;
 442     }
 443   }
 444
 445   Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
 446   if (diff > 0.1)
 447     amp = ymMaxY[0] ;
 448   //printf("Yves   : Amp %f, time %g\n",amp, time);
 449   //END YS
 450   return;
 451 }
 452
 453
 454
 455 //__________________________________________________________________
 456 Double_t
 457 AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t *x, Double_t *par)
 458 {
 459   // Matches version used in 2007 beam test
 460   //
 461   // Shape of the electronics raw reponse:
 462   // It is a semi-gaussian, 2nd order Gamma function of the general form
 463   //
 464   // xx = (t - t0 + tau) / tau  [xx is just a convenient help variable]
 465   // F = A * (xx**N * exp( N * ( 1 - xx) )   for xx >= 0
 466   // F = 0                                   for xx < 0
 467   //
 468   // parameters:
 469   // A:   par[0]   // Amplitude = peak value
 470   // t0:  par[1]
 471   // tau: par[2]
 472   // N:   par[3]
 473   // ped: par[4]
 474   //
 475   Double_t signal = 0.;
 476   Double_t tau    = par[2];
 477   Double_t n      = par[3];
 478   Double_t ped    = par[4];
 479   Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
 480
 481   if (xx <= 0)
 482     signal = ped ;
 483   else {
 484     signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
 485   }
 486   return signal ;
 487 }
 488