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 decodes the stream of ALTRO samples to extract
19 // the PHOS "digits" of current event. Uses fitting procedure
20 // to separate reasonable samples
23 // AliRawReader* rf = new AliRawReaderDate("2006run2211.raw");
24 // AliPHOSRawDecoder dc(rf);
25 // while (rf->NextEvent()) {
26 // dc.SetOldRCUFormat(kTRUE);
27 // dc.SubtractPedestals(kTRUE);
28 // while ( dc.NextDigit() ) {
29 // Int_t module = dc.GetModule();
30 // Int_t column = dc.GetColumn();
31 // Int_t row = dc.GetRow();
32 // Double_t amplitude = dc.GetEnergy();
33 // Double_t time = dc.GetTime();
34 // Bool_t IsLowGain = dc.IsLowGain();
39 // Author: Dmitri Peressounko
41 // --- ROOT system ---
52 // --- AliRoot header files ---
54 #include "AliPHOSRawDecoderv1.h"
55 #include "AliPHOSPulseGenerator.h"
58 ClassImp(AliPHOSRawDecoderv1)
60 //-----------------------------------------------------------------------------
61 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder(),
62 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
64 //Default constructor.
67 //-----------------------------------------------------------------------------
68 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader, AliAltroMapping **mapping):
69 AliPHOSRawDecoder(rawReader,mapping),
70 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
72 //Construct a decoder object.
73 //Is is user responsibility to provide next raw event
74 //using AliRawReader::NextEvent().
77 gMinuit = new TMinuit(100);
78 fSampleParamsHigh =new TArrayD(5) ;
79 fSampleParamsHigh->AddAt(4.25,0) ;
80 fSampleParamsHigh->AddAt(0.094,1) ;
81 fSampleParamsHigh->AddAt(0.0151,2) ;
82 fSampleParamsHigh->AddAt(0.0384,3) ;
83 fSampleParamsLow=new TArrayD(5) ;
84 fSampleParamsLow->AddAt(5.14,0) ;
85 fSampleParamsLow->AddAt(0.0970,1) ;
86 fSampleParamsLow->AddAt(0.0088,2) ;
87 fSampleParamsLow->AddAt(0.0346,3) ;
88 fToFit = new TList() ;
91 //-----------------------------------------------------------------------------
92 AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1()
96 delete fSampleParamsLow ;
99 if(fSampleParamsHigh){
100 delete fSampleParamsHigh ;
109 //-----------------------------------------------------------------------------
110 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ):
111 AliPHOSRawDecoder(phosDecoder),
112 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
115 fToFit = new TList() ;
116 fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
117 fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
120 //-----------------------------------------------------------------------------
121 AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecoder)
123 //Assignment operator.
125 // if(this != &phosDecoder) {
127 fToFit = new TList() ;
128 if(fSampleParamsLow){
129 fSampleParamsLow = phosDecoder.fSampleParamsLow ;
130 fSampleParamsHigh= phosDecoder.fSampleParamsHigh ;
133 fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
134 fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
139 //-----------------------------------------------------------------------------
140 Bool_t AliPHOSRawDecoderv1::NextDigit()
142 //Extract an energy deposited in the crystal,
143 //crystal' position (module,column,row),
144 //time and gain (high or low).
145 //First collects sample, then evaluates it and if it has
146 //reasonable shape, fits it with Gamma2 function and extracts
149 // TCanvas * c = (TCanvas *)gROOT->FindObjectAny("canvMy") ;
150 // TH1S * h = new TH1S("s","",200,0.5,200.5) ;
151 // TF1 * fff = new TF1("fff","[0]+[1]*((x-[2])+[3]*(x-[2])*(x-[2]))*(exp(-(x-[2])*[4])+[5]*exp(-(x-[2])*[6]))",0.,1000.) ;
153 AliCaloRawStream* in = fCaloStream;
160 Float_t baseLine = 1.0;
161 const Float_t nPreSamples = 10;
163 while ( in->Next() ) {
166 tLength = in->GetTimeLength();
167 if(tLength!=fSamples->GetSize()) {
169 fSamples = new TArrayI(tLength);
172 for(Int_t i=0; i<fSamples->GetSize(); i++){
173 fSamples->AddAt(0,i) ;
178 // Fit the full sample
179 if(in->IsNewHWAddress() && iBin>0) {
181 Double_t pedestal =0. ;
184 pedestal = (Double_t)(pedMean/nPed);
190 //first estimate if this sample looks like gamma2 function
195 for(Int_t i=0; i<fSamples->GetSize(); i++){
196 if(fSamples->At(i)>0){
197 Double_t de=fSamples->At(i)-pedestal ;
201 if(de>2 && tStart==0)
203 if(fSamples->At(i)>fEnergy)
204 fEnergy=fSamples->At(i) ;
209 aRMS=aRMS/cnts-aMean*aMean;
212 //IF sample has reasonable mean and RMS, try to fit it with gamma2
213 if(fEnergy>2.&& cnts >20 && aMean>0. && aRMS>2.){ //more or less reasonable sample
215 gMinuit->mncler(); // Reset Minuit's list of paramters
216 gMinuit->SetPrintLevel(-1) ; // No Printout
217 gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ;
218 // To set the address of the minimization function
222 fSampleParamsLow->AddAt(pedestal,4) ;
223 fToFit->AddFirst((TObject*)fSampleParamsLow) ;
226 fSampleParamsHigh->AddAt(pedestal,4) ;
227 fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
229 fToFit->AddLast((TObject*)fSamples) ;
231 gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare
233 gMinuit->mnparm(0, "t0", 1.*tStart, 0.1, 0, 0, ierflg) ;
235 // AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
238 return kTRUE ; //will scan further
240 Double_t amp0=(fEnergy-pedestal)*0.0032;
242 gMinuit->mnparm(1, "Energy", amp0 , 0.001*amp0, 0, 0, ierflg) ;
244 // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
247 return kTRUE ; //will scan further
250 Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
254 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
255 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
256 // gMinuit->SetMaxIterations(100);
257 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
259 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
263 gMinuit->GetParameter(0,t0, t0err) ;
264 gMinuit->GetParameter(1,efit, err) ;
268 a=fSampleParamsLow->At(0) ;
269 alpha=fSampleParamsLow->At(1) ;
272 a=fSampleParamsHigh->At(0) ;
273 alpha=fSampleParamsHigh->At(1) ;
279 // fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsLow->At(2),fSampleParamsLow->At(3)) ;
282 // fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsHigh->At(2),fSampleParamsHigh->At(3)) ;
284 // fff->Draw("same") ;
287 efit*=(2.*a+TMath::Sqrt(4.*a*a+alpha*alpha))/alpha/alpha*TMath::Exp(-1.+(alpha-TMath::Sqrt(4.*a*a+alpha*alpha))/2./a) ;
288 //printf("efit=%f, t0=%e +- %e, ped=%f \n",efit,t0,t0err,pedestal) ;
289 Double_t fmin,fedm,errdef ;
290 Int_t npari,nparx,istat;
291 gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
294 // printf("LowGain \n") ;
296 // printf("highGain \n") ;
298 //printf("fmin=%e \n",fmin) ;
301 if(1){ //fmin < 3.+0.3*efit ){ //Chi^2 of a good sample
308 fEnergy=0 ; //bad sample
312 fEnergy *= fPulseGenerator->GetRawFormatHighLowGainFactor(); // *16
314 fTime*=fPulseGenerator->GetRawFormatTimeTrigger() ;
316 if (fEnergy < baseLine) fEnergy = 0;
326 fLowGainFlag = in->IsLowGain();
327 fTime = fPulseGenerator->GetRawFormatTimeTrigger() * in->GetTime();
328 fModule = in->GetModule()+1;
329 fRow = in->GetRow() +1;
330 fColumn = in->GetColumn()+1;
333 // Fill array with samples
335 if(tLength-iBin < nPreSamples) {
336 pedMean += in->GetSignal();
339 fSamples->AddAt(in->GetSignal(),tLength-iBin);
340 // h->SetBinContent(tLength-iBin+1,in->GetSignal()) ;
346 //_____________________________________________________________________________
347 void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
349 // Calculates the Chi square for the samples minimization
350 // Number of parameters, Gradient, Chi squared, parameters, what to do
352 TList * toFit= (TList*)gMinuit->GetObjectFit() ;
353 TArrayD * params=(TArrayD*)toFit->At(0) ;
354 TArrayI * samples = (TArrayI*)toFit->At(1) ;
358 for(Int_t iparam = 0 ; iparam < 2 ; iparam++)
359 Grad[iparam] = 0 ; // Will evaluate gradient
361 Int_t nSamples=samples->GetSize() ; //Math::Min(70,samples->GetSize()) ;
364 Double_t a=params->At(0) ;
365 Double_t alpha=params->At(1) ;
366 Double_t b=params->At(2) ;
367 Double_t beta=params->At(3) ;
369 for(Int_t i = 0 ; i < nSamples ; i++) {
370 if(samples->At(i)==0 || samples->At(i)==1023) //zero or overflow
372 Double_t dt=i*1.-t0 ;
373 Double_t diff=float(samples->At(i))-Gamma2(dt,en,params) ;
374 Double_t w=0.1+0.005*i ; //Mean Pedestal RMS + rising modulation
377 if(iflag == 2){ // calculate gradient
379 Grad[0] += -2.*en*diff*((alpha*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-alpha*dt)+
380 b*(beta*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-beta*dt)) /w ; //derivative over t0
381 Grad[1] += -2.*diff*(dt+a*dt*dt)*(TMath::Exp(-alpha*dt)+
382 b*TMath::Exp(-dt*beta))/w ;
390 for(Int_t iparam = 0 ; iparam < 2 ; iparam++)
391 Grad[iparam] /= nSamples ;
396 //-----------------------------------------------------------------------------
397 Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t en,TArrayD * params){ //Function for fitting samples
399 //dt-time after start
401 //function parameters
403 Double_t ped=params->At(4) ;
405 return ped ; //pedestal
407 return ped+en*(dt+params->At(0)*dt*dt)*(TMath::Exp(-dt*params->At(1))+params->At(2)*TMath::Exp(-dt*params->At(3))) ;