1 /**************************************************************************
2 * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // This class extracts the signal parameters (energy, time, quality)
19 // from ALTRO samples. Energy is in ADC counts, time is in time bin units.
20 // Class uses FastFitting algorithm to fit sample and extract time and Amplitude
21 // and evaluate sample quality = (chi^2/NDF)/some parameterization providing
22 // efficiency close to 100%
25 // AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
26 // fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
27 // fitter->SetCalibData(fgCalibData) ;
28 // fitter->Eval(sig,sigStart,sigLength);
29 // Double_t amplitude = fitter.GetEnergy();
30 // Double_t time = fitter.GetTime();
31 // Bool_t isLowGain = fitter.GetCaloFlag()==0;
33 // Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
35 // --- ROOT system ---
45 // --- AliRoot header files ---
47 #include "AliPHOSCalibData.h"
48 #include "AliPHOSRawFitterv3.h"
49 #include "AliPHOSPulseGenerator.h"
51 ClassImp(AliPHOSRawFitterv3)
53 //-----------------------------------------------------------------------------
54 AliPHOSRawFitterv3::AliPHOSRawFitterv3():
57 //Default constructor.
60 //-----------------------------------------------------------------------------
61 AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
66 //-----------------------------------------------------------------------------
67 AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
68 AliPHOSRawFitterv0(phosFitter)
73 //-----------------------------------------------------------------------------
74 AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
76 //Assignment operator.
80 //-----------------------------------------------------------------------------
81 Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
83 //Extract an energy deposited in the crystal,
84 //crystal' position (module,column,row),
85 //time and gain (high or low).
86 //First collects sample, then evaluates it and if it has
87 //reasonable shape, fits it with Gamma2 function and extracts
90 if (fCaloFlag == 2 || fNBunches > 1) {
94 if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
103 const Float_t kBaseLine = 1.0;
104 const Int_t kPreSamples = 10;
107 //We tryed individual taus for each channel, but
108 //this approach seems to be unstable. Much better results are obtaned with
109 //fixed decay time for all channels.
110 const Double_t tau=22.18 ;
112 TArrayD samples(sigLength); // array of sample amplitudes
113 TArrayD times(sigLength); // array of sample time stamps
114 for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
116 pedMean += signal[i];
117 pedRMS += signal[i]*signal[i] ;
122 const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
123 Double_t pedestal = 0;
127 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
128 if(fPedestalRMS > 0.)
129 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
130 pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
131 fEnergy -= pedestal; // pedestal subtraction
137 //take pedestals from DB
138 pedestal = (Double_t) fAmpOffset ;
140 Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
141 Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
142 pedestal += truePed - altroSettings ;
145 // AliWarning(Form("Can not read data from OCDB")) ;
150 if (fEnergy < kBaseLine) fEnergy = 0;
152 fTime = sigStart-sigLength-3;
153 for (Int_t i=0; i<sigLength; i++) {
154 samples.AddAt(signal[i]-pedestal,sigLength-i-1);
155 times.AddAt(i/tau ,i);
158 //calculate time and energy
161 Int_t nMax = 0 ; //number of points in plato
165 for (Int_t i=0; i<sigLength; i++){
166 if(signal[i] > pedestal){
167 Double_t de = signal[i] - pedestal ;
173 if(signal[i] > maxAmp){
178 if(signal[i] == maxAmp){
184 if (maxBin==sigLength-1){//bad "rising" sample
191 fEnergy=Double_t(maxAmp)-pedestal ;
192 fOverflow =0 ; //look for plato on the top of sample
193 if (fEnergy>500 && //this is not fluctuation of soft sample
194 nMax>2){ //and there is a plato
200 aRMS = aRMS/wts - aMean*aMean;
203 //do not take too small energies
204 if (fEnergy < kBaseLine)
207 //do not test quality of too soft samples
208 if (fEnergy < maxEtoFit){
209 if (aRMS < 2.) //sigle peak
213 //Evaluate time of signal arriving
217 // if sample has reasonable mean and RMS, try to fit it with gamma2
218 //This method can not analyse overflow samples
224 Double_t a=0,b=0,c=0 ;
228 for(Int_t i=minI; i<sigLength; i++){
229 if(samples.At(i)<=0.)
231 Double_t t= times.At(i) ;
232 Double_t f02 = TMath::Exp(-2.*t);
233 Double_t f12 = t*f02;
234 Double_t f22 = t*f12;
236 Double_t f02d = -2.*f02;
237 Double_t f12d = f02 - 2.*f12;
238 Double_t f22d = 2.*(f12 - f22);
239 a += f02d * samples.At(i) ;
240 b -= f12d * samples.At(i) ;
241 c += f22d * samples.At(i) ;
246 if(b==0.){ //no roots
248 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
249 printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
255 Double_t amp=0.,den=0.; ;
256 for(Int_t i=minI; i<sigLength; i++){
257 if(samples.At(i)<=0.)
259 Double_t dt = times.At(i) - t1;
260 Double_t f = dt*dt*TMath::Exp(-2.*dt);
261 amp += f*samples.At(i);
264 if(den>0.0) amp /= den;
267 for(Int_t i=minI; i<sigLength; i++){
268 if(samples.At(i)<=0.)
270 Double_t t = times.At(i)- t1;
271 Double_t dy = samples.At(i)- amp*t*t*TMath::Exp(-2.*t) ;
275 fEnergy = amp*TMath::Exp(-2.);
276 fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
279 Double_t det = b*b - a*c;
280 if(det>=1.e-6 && det<0.0) {
281 det = 0.0; //rounding error
283 if(det<0.){ //Problem
285 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
286 printf(" det=%e \n",det) ;
292 det = TMath::Sqrt(det);
293 Double_t t1 = (-b + det) / a;
294 // Double_t t2 = (-b - det) / a; //second root is wrong one
295 Double_t amp1=0., den1=0. ;
296 for(Int_t i=minI; i<sigLength; i++){
297 if(samples.At(i)<=0.)
299 Double_t dt1 = times.At(i) - t1;
300 Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
301 amp1 += f01*samples.At(i);
304 if(den1>0.0) amp1 /= den1;
305 Double_t chi1=0.; // chi2=0. ;
306 for(Int_t i=minI; i<sigLength; i++){
307 if(samples.At(i)<=0.)
309 Double_t dt1 = times.At(i)- t1;
310 Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
313 fEnergy=amp1*TMath::Exp(-2.) ; ;
315 fQuality=chi1/sigLength ;
318 //Impose cut on quality
319 fQuality/=2.+0.004*fEnergy*fEnergy ;
321 //Draw corrupted samples
322 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
324 if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
325 // if(!fOverflow ){ //Draw only bad samples
326 printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
327 TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
328 if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
330 for (Int_t i=0; i<sigLength; i++) {
331 h->SetBinContent(i+1,samples.At(i)+pedestal) ;
333 TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
334 fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
335 fffit->SetLineColor(2) ;
336 TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
338 can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
339 can->SetFillColor(0) ;
340 can->SetFillStyle(0) ;
342 can->SetBorderSize(0);
346 TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
349 spectrum_1->Range(0,0,1,1);
350 spectrum_1->SetFillColor(0);
351 spectrum_1->SetFillStyle(4000);
352 spectrum_1->SetBorderSize(1);
353 spectrum_1->SetBottomMargin(0.012);
354 spectrum_1->SetTopMargin(0.03);
355 spectrum_1->SetLeftMargin(0.10);
356 spectrum_1->SetRightMargin(0.05);
359 snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
362 fffit->Draw("same") ;
364 snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
365 TFile fout("samples_bad.root","update") ;
370 TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
371 spectrum_2->SetFillColor(0) ;
372 spectrum_2->SetFillStyle(0) ;
373 spectrum_2->SetGridy() ;
375 spectrum_2->Range(0,0,1,1);
376 spectrum_2->SetFillColor(0);
377 spectrum_2->SetBorderSize(1);
378 spectrum_2->SetTopMargin(0.01);
379 spectrum_2->SetBottomMargin(0.25);
380 spectrum_2->SetLeftMargin(0.10);
381 spectrum_2->SetRightMargin(0.05);
384 TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
385 if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
387 for (Int_t i=0; i<sigLength; i++) {
388 hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
393 printf("Press <enter> to continue\n") ;