//We tryed individual taus for each channel, but
//this approach seems to be unstable. Much better results are obtaned with
//fixed decay time for all channels.
- const Double_t tau=21. ;
+ const Double_t tau=22.18 ;
- TArrayD *samples = new TArrayD(sigLength); // array of sample amplitudes
- TArrayD *times = new TArrayD(sigLength); // array of sample time stamps
- for (Int_t i=0; i<sigLength; i++) {
- if (i<kPreSamples) {
- nPed++;
- pedMean += signal[i];
- pedRMS += signal[i]*signal[i] ;
- }
- samples->AddAt(signal[i],sigLength-i-1);
- times ->AddAt(i/tau ,i);
+ TArrayD samples(sigLength); // array of sample amplitudes
+ TArrayD times(sigLength); // array of sample time stamps
+ for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
+ nPed++;
+ pedMean += signal[i];
+ pedRMS += signal[i]*signal[i] ;
}
fEnergy = -111;
fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
if(fPedestalRMS > 0.)
fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
- fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
+ pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
+ fEnergy -= pedestal; // pedestal subtraction
}
else
return kFALSE;
pedestal += truePed - altroSettings ;
}
else{
- AliWarning(Form("Can not read data from OCDB")) ;
+// AliWarning(Form("Can not read data from OCDB")) ;
}
fEnergy-=pedestal ;
}
if (fEnergy < kBaseLine) fEnergy = 0;
//Evaluate time
- Int_t iStart = 0;
- while(iStart<sigLength && samples->At(iStart)-pedestal <kBaseLine) iStart++ ;
- fTime = sigStart+iStart;
+ fTime = sigStart-sigLength-3;
+ for (Int_t i=0; i<sigLength; i++) {
+ samples.AddAt(signal[i]-pedestal,sigLength-i-1);
+ times.AddAt(i/tau ,i);
+ }
//calculate time and energy
Int_t maxBin=0 ;
Double_t aMean =0. ;
Double_t aRMS =0. ;
Double_t wts =0 ;
- Int_t tStart=0 ;
-
for (Int_t i=0; i<sigLength; i++){
if(signal[i] > pedestal){
Double_t de = signal[i] - pedestal ;
aRMS += de*i*i ;
wts += de;
}
- if(de > 2 && tStart==0)
- tStart = i ;
if(signal[i] > maxAmp){
maxAmp = signal[i];
nMax=0;
maxBin = i ;
}
- if(signal[i] == maxAmp)
+ if(signal[i] == maxAmp){
nMax++;
+ }
}
}
//do not test quality of too soft samples
if (fEnergy < maxEtoFit){
- fTime = tStart;
if (aRMS < 2.) //sigle peak
fQuality = 999. ;
else
fQuality = 0. ;
+ //Evaluate time of signal arriving
return kTRUE ;
}
// if sample has reasonable mean and RMS, try to fit it with gamma2
+ //This method can not analyse overflow samples
+ if(fOverflow){
+ fQuality = 99. ;
+ return kTRUE ;
+ }
// First estimate t0
Double_t a=0,b=0,c=0 ;
- for(Int_t i=0; i<sigLength; i++){
- if(samples->At(i)<pedestal)
+ Int_t minI=0 ;
+ if (fPedSubtract)
+ minI=kPreSamples ;
+ for(Int_t i=minI; i<sigLength; i++){
+ if(samples.At(i)<=0.)
continue ;
- Double_t t= times->At(i) ;
+ Double_t t= times.At(i) ;
Double_t f02 = TMath::Exp(-2.*t);
Double_t f12 = t*f02;
Double_t f22 = t*f12;
Double_t f02d = -2.*f02;
Double_t f12d = f02 - 2.*f12;
Double_t f22d = 2.*(f12 - f22);
- a += f02d * (samples->At(i)-pedestal) ;
- b -= f12d * (samples->At(i)-pedestal) ;
- c += f22d * (samples->At(i)-pedestal) ;
+ a += f02d * samples.At(i) ;
+ b -= f12d * samples.At(i) ;
+ c += f22d * samples.At(i) ;
}
- //Find 2 roots
- Double_t det = b*b - a*c;
- if(det>=1.e-6 && det<0.0) {
- det = 0.0; //rounding error
- }
- if(det<0.){ //Problem
- fQuality = 1500. ;
- if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
- printf(" det=%e \n",det) ;
- goto plot ;
+ //Find roots
+ if(a==0.){
+ if(b==0.){ //no roots
+ fQuality = 2000. ;
+ if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
+ printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
+ goto plot ;
+ }
+ return kTRUE ;
}
- return kTRUE ;
+ Double_t t1=-c/b ;
+ Double_t amp=0.,den=0.; ;
+ for(Int_t i=minI; i<sigLength; i++){
+ if(samples.At(i)<=0.)
+ continue ;
+ Double_t dt = times.At(i) - t1;
+ Double_t f = dt*dt*TMath::Exp(-2.*dt);
+ amp += f*samples.At(i);
+ den += f*f;
+ }
+ if(den>0.0) amp /= den;
+ // chi2 calculation
+ fQuality=0.;
+ for(Int_t i=minI; i<sigLength; i++){
+ if(samples.At(i)<=0.)
+ continue ;
+ Double_t t = times.At(i)- t1;
+ Double_t dy = samples.At(i)- amp*t*t*TMath::Exp(-2.*t) ;
+ fQuality += dy*dy;
+ }
+ fTime+=t1*tau ;
+ fEnergy = amp*TMath::Exp(-2.);
+ fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
}
+ else{
+ Double_t det = b*b - a*c;
+ if(det>=1.e-6 && det<0.0) {
+ det = 0.0; //rounding error
+ }
+ if(det<0.){ //Problem
+ fQuality = 1500. ;
+ if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
+ printf(" det=%e \n",det) ;
+ goto plot ;
+ }
+ return kTRUE ;
+ }
- if(det>0.0 && a!=0.) {
det = TMath::Sqrt(det);
Double_t t1 = (-b + det) / a;
// Double_t t2 = (-b - det) / a; //second root is wrong one
Double_t amp1=0., den1=0. ;
- for(Int_t i=0; i<sigLength; i++){
- Double_t dt1 = times->At(i) - t1;
- Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
- amp1 += f01*(samples->At(i)-pedestal);
- den1 += f01*f01;
- }
- if(den1>0.0) amp1 /= den1;
- Double_t chi1=0.; // chi2=0. ;
- for(Int_t i=0; i<sigLength; i++){
- Double_t dt1 = times->At(i)- t1;
- Double_t dy1 = samples->At(i)-pedestal- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
- chi1 += dy1*dy1;
- }
- fEnergy=amp1*TMath::Exp(-2.) ; ;
- fTime=t1*tau ;
- fQuality=chi1/sigLength ;
+ for(Int_t i=minI; i<sigLength; i++){
+ if(samples.At(i)<=0.)
+ continue ;
+ Double_t dt1 = times.At(i) - t1;
+ Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
+ amp1 += f01*samples.At(i);
+ den1 += f01*f01;
+ }
+ if(den1>0.0) amp1 /= den1;
+ Double_t chi1=0.; // chi2=0. ;
+ for(Int_t i=minI; i<sigLength; i++){
+ if(samples.At(i)<=0.)
+ continue ;
+ Double_t dt1 = times.At(i)- t1;
+ Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
+ chi1 += dy1*dy1;
+ }
+ fEnergy=amp1*TMath::Exp(-2.) ; ;
+ fTime+=t1*tau ;
+ fQuality=chi1/sigLength ;
}
- else {
- Double_t t1 ;
- if(a!=0.)
- t1=-b/a ;
- else
- t1=-c/b ;
- Double_t amp=0.,den=0.; ;
- for(Int_t i=0; i<sigLength; i++){
- Double_t dt = times->At(i) - t1;
- Double_t f = dt*dt*TMath::Exp(-2.*dt);
- amp += f*samples->At(i);
- den += f*f;
- }
- if(den>0.0) amp /= den;
- // chi2 calculation
- fQuality=0.;
- for(Int_t i=0; i<sigLength; i++){
- Double_t t = times->At(i)- t1;
- Double_t dy = samples->At(i)- amp*t*t*TMath::Exp(-2.*t) ;
- fQuality += dy*dy;
- }
- fTime=t1*tau ;
- fEnergy = amp*TMath::Exp(-2.);
- fQuality/= sigLength ;
- }
//Impose cut on quality
fQuality/=2.+0.004*fEnergy*fEnergy ;
//Draw corrupted samples
if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
plot:
- if(fQuality >1. && !fOverflow ){ //Draw only bad samples
+ if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
+// if(!fOverflow ){ //Draw only bad samples
printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
- if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
+ if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
h->Reset() ;
for (Int_t i=0; i<sigLength; i++) {
- h->SetBinContent(i,samples->At(i)) ;
+ h->SetBinContent(i+1,samples.At(i)+pedestal) ;
}
TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
fffit->SetLineColor(2) ;
- TCanvas * c = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
- if(!c){
- c = new TCanvas("cSamples","cSamples",10,10,400,400) ;
- c->SetFillColor(0) ;
- c->SetFillStyle(0) ;
- c->Range(0,0,1,1);
- c->SetBorderSize(0);
+ TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
+ if(!can){
+ can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
+ can->SetFillColor(0) ;
+ can->SetFillStyle(0) ;
+ can->Range(0,0,1,1);
+ can->SetBorderSize(0);
}
- c->cd() ;
+ can->cd() ;
TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
spectrum_1->Draw();
spectrum_1->SetRightMargin(0.05);
char title[155] ;
- sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
+ snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
h->SetTitle(title) ;
h->Draw() ;
fffit->Draw("same") ;
- sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
+ snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
TFile fout("samples_bad.root","update") ;
h->Write(title);
fout.Close() ;
- c->cd() ;
+ can->cd() ;
TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
spectrum_2->SetFillColor(0) ;
spectrum_2->SetFillStyle(0) ;
spectrum_2->cd() ;
TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
- if(!hd) hd = new TH1I("hd","Samples",50,0.,50.) ;
+ if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
hd->Reset() ;
for (Int_t i=0; i<sigLength; i++) {
- hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
+ hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
}
hd->Draw();
-
- c->Update() ;
+/*
+ can->Update() ;
printf("Press <enter> to continue\n") ;
getchar();
-
+*/
delete fffit ;
delete spectrum_1 ;
}
}
- delete samples ;
- delete times ;
return kTRUE;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff