PHOS/AliPHOSRawFitterv4.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 // This class extracts the signal parameters (energy, time, quality)
  18 // from ALTRO samples. Energy is in ADC counts, time is in time bin units.
  19 // If sample is not in saturation, a coarse algorithm is used (a-la AliPHOSRawFitterv0)
  20 // If sample is in saturation, the unsaturated part of the sample is fit a-la AliPHOSRawFitterv1
  21 //
  22 //     AliPHOSRawFitterv4 *fitterv4=new AliPHOSRawFitterv4();
  23 //     fitterv4->SetChannelGeo(module,cellX,cellZ,caloFlag);
  24 //     fitterv4->SetCalibData(fgCalibData) ;
  25 //     fitterv4->Eval(sig,sigStart,sigLength);
  26 //     Double_t amplitude = fitterv4.GetEnergy();
  27 //     Double_t time      = fitterv4.GetTime();
  28 //     Bool_t   isLowGain = fitterv4.GetCaloFlag()==0;
  29
  30 // Author: Dmitri Peressounko
  31
  32 // --- ROOT system ---
  33 #include "TArrayI.h"
  34 #include "TMath.h"
  35 #include "TObject.h"
  36 #include "TArrayD.h"
  37 #include "TList.h"
  38 #include "TMinuit.h"
  39
  40
  41 //Used for debug
  42 //#include "TROOT.h"
  43 //#include "TH1.h"
  44 //#include "TCanvas.h"
  45 //#include "TPad.h"
  46 //#include "TF1.h"
  47
  48
  49
  50 // --- AliRoot header files ---
  51 #include "AliPHOSRawFitterv4.h"
  52 #include "AliPHOSCalibData.h"
  53 #include "AliLog.h"
  54
  55 ClassImp(AliPHOSRawFitterv4)
  56
  57 //-----------------------------------------------------------------------------
  58 AliPHOSRawFitterv4::AliPHOSRawFitterv4():
  59   AliPHOSRawFitterv1(),
  60   fFitHighGain(0)
  61 {
  62   //Default constructor
  63 }
  64
  65 //-----------------------------------------------------------------------------
  66 AliPHOSRawFitterv4::~AliPHOSRawFitterv4()
  67 {
  68 }
  69
  70 //-----------------------------------------------------------------------------
  71 AliPHOSRawFitterv4::AliPHOSRawFitterv4(const AliPHOSRawFitterv4 &phosFitter ):
  72   AliPHOSRawFitterv1((AliPHOSRawFitterv1)phosFitter),
  73   fFitHighGain(0)
  74 {
  75   //Copy constructor
  76   fFitHighGain=phosFitter.fFitHighGain;
  77 }
  78
  79
  80 //-----------------------------------------------------------------------------
  81
  82 Bool_t AliPHOSRawFitterv4::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
  83 {
  84   // Calculate signal parameters (energy, time, quality) from array of samples
  85   // Energy is a maximum sample minus pedestal 9
  86   // Time is the first time bin
  87   // Signal overflows is there are at least 3 samples of the same amplitude above 900
  88
  89   fOverflow= kFALSE ;
  90   fEnergy  = 0;
  91   if (fNBunches > 1) {
  92     fQuality = 1000;
  93     return kTRUE;
  94   }
  95
  96   const Float_t kBaseLine   = 1.0;
  97   const Int_t   kPreSamples = 10;
  98
  99   Float_t  pedMean   = 0;
 100   Float_t  pedRMS    = 0;
 101   Int_t    nPed      = 0;
 102   UShort_t maxSample = 0;
 103   Int_t    nMax      = 0;
 104
 105   for (Int_t i=0; i<sigLength; i++) {
 106     if (i>sigLength-kPreSamples) { //inverse signal time order
 107       nPed++;
 108       pedMean += signal[i];
 109       pedRMS  += signal[i]*signal[i] ;
 110     }
 111     if(signal[i] > maxSample ){ maxSample = signal[i]; nMax=0;}
 112     if(signal[i] == maxSample) nMax++;
 113
 114   }
 115
 116   fEnergy = (Double_t)maxSample;
 117   if (maxSample > 850 && nMax > 2) fOverflow = kTRUE;
 118
 119   Double_t pedestal = 0 ;
 120   if (fPedSubtract) {
 121     if (nPed > 0) {
 122       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 123       if(fPedestalRMS > 0.)
 124         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 125       pedestal = (Double_t)(pedMean/nPed);
 126     }
 127     else
 128       return kFALSE;
 129   }
 130   else {
 131     //take pedestals from DB
 132     pedestal = (Double_t) fAmpOffset ;
 133     if (fCalibData) {
 134       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
 135       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
 136       pedestal += truePed - altroSettings ;
 137     }
 138     else{
 139       AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
 140     }
 141   }
 142   fEnergy-=pedestal ;
 143   if (fEnergy < kBaseLine) fEnergy = 0;
 144
 145
 146   if(fOverflow && ((fCaloFlag == 0) || fFitHighGain)){ //by default fit LowGain only
 147     TArrayI *samples = new TArrayI(sigLength); // array of sample amplitudes
 148     TArrayI *times   = new TArrayI(sigLength); // array of sample time stamps
 149     Bool_t result = kTRUE ;
 150     for (Int_t i=0; i<sigLength; i++) {
 151       samples->AddAt(signal[i]-pedestal,sigLength-i-1);
 152       times  ->AddAt(i ,i);
 153     }
 154     //Prepare fit parameters
 155     //Evaluate time
 156     Int_t iStart = 0;
 157     while(iStart<sigLength && samples->At(iStart) <kBaseLine) iStart++ ;
 158     if      (fCaloFlag == 0){ // Low gain
 159       fSampleParamsLow->AddAt(pedestal,4) ;
 160       fSampleParamsLow->AddAt(double(maxSample),5) ;
 161       fSampleParamsLow->AddAt(double(iStart),6) ;
 162       fToFit->AddFirst((TObject*)fSampleParamsLow) ;
 163     }
 164     else if (fCaloFlag == 1){ // High gain
 165       fSampleParamsHigh->AddAt(pedestal,4) ;
 166       fSampleParamsHigh->AddAt(double(maxSample),5) ;
 167       fSampleParamsHigh->AddAt(double(iStart),6) ;
 168       fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
 169     }
 170     result=EvalWithFitting(samples,times);
 171     delete samples ;
 172     delete times ;
 173
 174     return result;
 175
 176   }
 177
 178
 179   //Sample withour overflow - Evaluate time
 180   fTime = sigStart-sigLength-3;
 181   const Int_t nLine= 6 ;        //Parameters of fitting
 182   const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
 183   const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
 184   // Avoid too low peak:
 185   if(fEnergy < eMinTOF){
 186      return kTRUE;
 187   }
 188   // Find index posK (kLevel is a level of "timestamp" point Tk):
 189   Int_t posK =sigLength-1 ; //last point before crossing k-level
 190   Double_t levelK = pedestal + kAmp*fEnergy;
 191   while(signal[posK] <= levelK && posK>=0){
 192      posK-- ;
 193   }
 194   posK++ ;
 195
 196   if(posK == 0 || posK==sigLength-1){
 197     return kTRUE;
 198   }
 199
 200   // Find crosing point by solving linear equation (least squares method)
 201   Int_t np = 0;
 202   Int_t iup=posK-1 ;
 203   Int_t idn=posK ;
 204   Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
 205   Double_t x,y ;
 206
 207   while(np<nLine){
 208     //point above crossing point
 209     if(iup>=0){
 210       x = sigLength-iup-1;
 211       y = signal[iup];
 212       sx += x;
 213       sy += y;
 214       sxx += (x*x);
 215       sxy += (x*y);
 216       np++ ;
 217       iup-- ;
 218     }
 219     //Point below crossing point
 220     if(idn<sigLength){
 221       if(signal[idn]<pedestal){
 222         idn=sigLength-1 ; //do not scan further
 223         idn++ ;
 224         continue ;
 225       }
 226       x = sigLength-idn-1;
 227       y = signal[idn];
 228       sx += x;
 229       sy += y;
 230       sxx += (x*x);
 231       sxy += (x*y);
 232       np++;
 233       idn++ ;
 234     }
 235     if(idn>=sigLength && iup<0){
 236       break ; //can not fit futher
 237     }
 238   }
 239
 240   Double_t det = np*sxx - sx*sx;
 241   if(det == 0){
 242     return kTRUE;
 243   }
 244   Double_t c1 = (np*sxy - sx*sy)/det;  //slope
 245   Double_t c0 = (sy-c1*sx)/np; //offset
 246   if(c1 == 0){
 247     return kTRUE;
 248   }
 249
 250   // Find where the line cross kLevel:
 251   fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
 252   return kTRUE;
 253
 254 }
 255 //===================================================
 256 Bool_t AliPHOSRawFitterv4::EvalWithFitting(TArrayI*samples, TArrayI * times){
 257
 258   // if sample has reasonable mean and RMS, try to fit it with gamma2
 259   const Float_t sampleMaxHG=102.332 ;
 260   const Float_t sampleMaxLG=277.196 ;
 261
 262   gMinuit->mncler();                     // Reset Minuit's list of paramters
 263   gMinuit->SetPrintLevel(-1) ;           // No Printout
 264   gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ;
 265   // To set the address of the minimization function
 266
 267   fToFit->Clear("nodelete") ;
 268   Double_t b=0,bmin=0,bmax=0 ;
 269   if      (fCaloFlag == 0){ // Low gain
 270     b=fSampleParamsLow->At(2) ;
 271     bmin=0.5 ;
 272     bmax=10. ;
 273     fToFit->AddFirst((TObject*)fSampleParamsLow) ;
 274   }
 275   else if (fCaloFlag == 1){ // High gain
 276     b=fSampleParamsHigh->At(2) ;
 277     bmin=0.05 ;
 278     bmax=0.4 ;
 279     fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
 280   }
 281   fToFit->AddLast((TObject*)samples) ;
 282   fToFit->AddLast((TObject*)times) ;
 283
 284
 285   gMinuit->SetObjectFit((TObject*)fToFit) ;         // To tranfer pointer to UnfoldingChiSquare
 286   Int_t ierflg ;
 287   gMinuit->mnparm(0, "t0",  1., 0.01, -50., 50., ierflg) ;
 288   if(ierflg != 0){
 289     //    AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
 290     fEnergy =   0. ;
 291     fTime   =-999. ;
 292     fQuality= 999. ;
 293     return kTRUE ; //will scan further
 294   }
 295   Double_t amp0=0;
 296   if      (fCaloFlag == 0) // Low gain
 297     amp0 = fEnergy/sampleMaxLG;
 298   else if (fCaloFlag == 1) // High gain
 299     amp0 = fEnergy/sampleMaxHG;
 300
 301   gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
 302   if(ierflg != 0){
 303     //    AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
 304     fEnergy =   0. ;
 305     fTime   =-999. ;
 306     fQuality= 999. ;
 307     return kTRUE ; //will scan further
 308   }
 309
 310   gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ;
 311   if(ierflg != 0){
 312     //        AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
 313     fEnergy =   0. ;
 314     fTime   =-999. ;
 315     fQuality= 999. ;
 316     return kTRUE ; //will scan further
 317   }
 318
 319   Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 320   //  depends on it.
 321   Double_t p1 = 1.0 ;
 322   Double_t p2 = 0.0 ;
 323   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 324   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 325   //////        gMinuit->SetMaxIterations(100);
 326   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 327
 328   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 329
 330   Double_t err=0.,t0err=0. ;
 331   Double_t t0=0.,efit=0. ;
 332   gMinuit->GetParameter(0,t0, t0err) ;
 333   gMinuit->GetParameter(1,efit, err) ;
 334
 335   Double_t bfit=0., berr=0. ;
 336   gMinuit->GetParameter(2,bfit,berr) ;
 337
 338    //Calculate total energy
 339   //this is parameterization of dependence of pulse height on parameter b
 340   if(fCaloFlag == 0) // Low gain
 341     fEnergy = efit*(99.54910 + 78.65038*bfit) ;
 342   else if(fCaloFlag == 1) // High gain
 343     fEnergy=efit*(80.33109 + 128.6433*bfit) ;
 344
 345   if(fEnergy < 0. || fEnergy > 10000.){
 346     //set energy to previously found max
 347     fTime   =-999.;
 348     fQuality= 999 ;
 349     fEnergy=0. ;
 350     return kTRUE;
 351   }
 352
 353   //evaluate fit quality
 354   Double_t fmin=0.,fedm=0.,errdef=0. ;
 355   Int_t npari,nparx,istat;
 356   gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
 357   fQuality = fmin/samples->GetSize() ;
 358   //compare quality with some parameterization
 359   if      (fCaloFlag == 0) // Low gain
 360     fQuality /= 2.00 + 0.0020*fEnergy ;
 361   else if (fCaloFlag == 1) // High gain
 362     fQuality /= 0.75 + 0.0025*fEnergy ;
 363
 364   fTime  += t0 - 4.024*bfit ;
 365
 366   if(fQuality==0.){//no points to fit)
 367     fTime   =-999.;
 368     fQuality= 1999 ;
 369     fEnergy=0. ;
 370   }
 371
 372 /*
 373   if(1){
 374     TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
 375     if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
 376     h->Reset() ;
 377     for (Int_t i=0; i<samples->GetSize(); i++) {
 378         h->SetBinContent(i+1,float(samples->At(i))) ;
 379     }
 380     TF1 * fffit = new TF1("fffit","[0]*(abs(x-[1])^[3]*exp(-[2]*(x-[1]))+[4]*(x-[1])*(x-[1])*exp(-[5]*(x-[1])))",0.,60.) ;
 381     TArrayD * params=(TArrayD*)fToFit->At(0) ;
 382     Double_t n=params->At(0) ;
 383     Double_t alpha=params->At(1) ;
 384     Double_t beta=params->At(3) ;
 385     fffit->SetParameters(efit,t0,alpha,n,bfit,beta) ;
 386     fffit->SetLineColor(2) ;
 387     TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
 388     if(!can){
 389         can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
 390         can->SetFillColor(0) ;
 391         can->SetFillStyle(0) ;
 392         can->Range(0,0,1,1);
 393         can->SetBorderSize(0);
 394     }
 395     can->cd() ;
 396
 397 //      TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
 398       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.001,0.999,0.999);
 399       spectrum_1->Draw();
 400       spectrum_1->cd();
 401       spectrum_1->Range(0,0,1,1);
 402       spectrum_1->SetFillColor(0);
 403       spectrum_1->SetFillStyle(4000);
 404       spectrum_1->SetBorderSize(1);
 405       spectrum_1->SetBottomMargin(0.012);
 406       spectrum_1->SetTopMargin(0.03);
 407       spectrum_1->SetLeftMargin(0.10);
 408       spectrum_1->SetRightMargin(0.05);
 409
 410       char title[155] ;
 411       snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
 412       h->SetTitle(title) ;
 413 //      h->Fit(fffit,"","",0.,51.) ;
 414       h->Draw() ;
 415       fffit->Draw("same") ;
 416
 417
 418 //       can->cd() ;
 419 //       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
 420 //       spectrum_2->SetFillColor(0) ;
 421 //       spectrum_2->SetFillStyle(0) ;
 422 //       spectrum_2->SetGridy() ;
 423 //       spectrum_2->Draw();
 424 //       spectrum_2->Range(0,0,1,1);
 425 //       spectrum_2->SetFillColor(0);
 426 //       spectrum_2->SetBorderSize(1);
 427 //       spectrum_2->SetTopMargin(0.01);
 428 //       spectrum_2->SetBottomMargin(0.25);
 429 //       spectrum_2->SetLeftMargin(0.10);
 430 //       spectrum_2->SetRightMargin(0.05);
 431 //       spectrum_2->cd() ;
 432 //
 433 //       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
 434 //       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
 435 //       hd->Reset() ;
 436 //       for (Int_t i=0; i<samples->GetSize(); i++) {
 437 //         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
 438 //       }
 439 //       hd->Draw();
 440
 441       can->Update() ;
 442       printf("Press <enter> to continue\n") ;
 443       getchar();
 444
 445
 446       delete fffit ;
 447       delete spectrum_1 ;
 448 //      delete spectrum_2 ;
 449     }
 450 */
 451
 452   return kTRUE;
 453
 454 }