Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, 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 : ///////////////////////////////////////////////////////////////////////////
18 : /// \file AliFFTsmoother.cxx
19 : /// \class AliFFTsmmother
20 : /// FFT smoothing algorithm. Transform input histogram removing outlier frequency or smoothing frequencies
21 : ///
22 :
23 :
24 : #include "TH1.h"
25 : #include "TGraph.h"
26 : #include "TTreeStream.h"
27 : #include "TVirtualFFT.h"
28 : #include "TMath.h"
29 : #include "TVector.h"
30 : #include "TStatToolkit.h"
31 : #include "AliFFTsmoother.h"
32 :
33 :
34 : TGraph * AliFFTsmoother::ReplaceOutlierFrequenciesMedian(TH1 *hinput, Double_t outlierCut, Int_t medianRange, Float_t smoothSigma, Int_t lowBand, TTreeSRedirector * pcstream)
35 : {
36 : ///
37 : /// FFT smoothing algorithm. Transform input histogram removing outlier frequency
38 : /// Code used for the analysis of the CE electron transparency measurement
39 : /// authors: Marian + Mesut
40 : /// Current algorithm is sensitive to frequency "aliasing"
41 : /// see detailed discussin in the UnitTest routines
42 : /// parasitic frequency and input histogram binning should be in "sync"
43 : /// e.g having repetetive structure 4 bins nx4 bins hsould be used
44 : ///
45 : /// Parameters:
46 : /// hinput - input histogram
47 : /// outlierCut - nsigma value abouve which the frequency is replaced |f_{i}-(f_{i+1)+f_{i-1})/2|<outlierCut*sigma
48 : /// lowBand - low freuency band (untoched)
49 : /// pcstream - ption to dump the intermedait results FFT into tree for further analysis
50 : ///
51 : /// Algorithm:
52 : /// 1.) Make FFT of the input histogram
53 : /// 2.) Replace outlier frequencies
54 : /// 3.) Make a back FFT transformation
55 : /// 4.) In case specified streamer dump intermediat results for the checking purposes
56 0 : TGraph grInput(hinput);
57 0 : Int_t fftLength = hinput->GetNbinsX();
58 : TH1D *hm = 0;
59 0 : TVirtualFFT::SetTransform(0);
60 0 : hm = (TH1D*)hinput->FFT(hm, "MAG");
61 0 : hm->SetTitle("Magnitude of the 1st transform");
62 :
63 : //
64 : // 1.) Make FFT of the input histogram
65 : //
66 0 : TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
67 0 : TVectorD reFull(fftLength);
68 0 : TVectorD imFull(fftLength);
69 0 : TVectorD magFull(fftLength);
70 0 : TVectorD magS(fftLength);
71 0 : TVectorD phaseFull(fftLength);
72 0 : fft->GetPointsComplex(reFull.GetMatrixArray(),imFull.GetMatrixArray());
73 0 : for (Int_t ipoint=1; ipoint<fftLength/2-1; ipoint++){
74 0 : magFull[ipoint]=TMath::Sqrt(reFull[ipoint]*reFull[ipoint]+imFull[ipoint]*imFull[ipoint]);
75 0 : phaseFull[ipoint]=TMath::ATan2(imFull[ipoint], reFull[ipoint]);
76 : }
77 : //
78 0 : TVectorD reMod(fftLength);
79 0 : TVectorD imMod(fftLength);
80 0 : TVectorD vecDiff(fftLength);
81 : //
82 : // 2.) identify and replace outlier frequencies
83 : //
84 0 : TVectorD magMed(fftLength);
85 0 : for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
86 0 : Int_t dedge=TMath::Min(ipoint, fftLength/2-ipoint-1);
87 0 : Int_t window= TMath::Min(dedge, medianRange);
88 0 : magMed[ipoint]=TMath::Median(1+2*window, &(magFull.GetMatrixArray()[ipoint-window]));
89 0 : vecDiff[ipoint]=0;
90 0 : if (magMed[ipoint]>0){
91 0 : vecDiff[ipoint]=(magFull[ipoint]-magMed[ipoint]);
92 0 : }
93 : }
94 0 : for (Int_t ipoint=1; ipoint<fftLength/2-1; ipoint++){
95 0 : vecDiff[ipoint]= (magFull[ipoint]-(magMed[ipoint-1]+magMed[ipoint+1])*0.5);
96 : }
97 : //
98 0 : Double_t meanT, rmsT;
99 0 : TStatToolkit::EvaluateUni(fftLength/2,vecDiff.GetMatrixArray(), meanT, rmsT, 0.95*(fftLength/2));
100 0 : if (smoothSigma<=0){
101 0 : for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
102 0 : reMod[ipoint]=reFull[ipoint];
103 0 : imMod[ipoint]=imFull[ipoint];
104 0 : if (ipoint<lowBand) continue;
105 0 : if (ipoint>=fftLength/2-lowBand) continue;
106 0 : if (TMath::Abs(vecDiff[ipoint]-meanT)>rmsT*outlierCut){
107 0 : reMod[ipoint]=magMed[ipoint]*TMath::Cos(phaseFull[ipoint]);
108 0 : imMod[ipoint]=magMed[ipoint]*TMath::Sin(phaseFull[ipoint]);
109 0 : }
110 : }
111 0 : }else{
112 0 : for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
113 0 : reMod[ipoint]=reFull[ipoint];
114 0 : imMod[ipoint]=imFull[ipoint];
115 0 : if (ipoint<lowBand) continue;
116 0 : if (ipoint>=fftLength/2-lowBand) continue;
117 : //
118 : Double_t sumM=0, sumW=0;
119 0 : for (Int_t dpoint=-4.*smoothSigma; dpoint<4.*smoothSigma; dpoint++){
120 0 : if (dpoint+ipoint<0) continue;
121 0 : if (dpoint+ipoint>fftLength/2-1) continue;
122 0 : if (TMath::Abs(vecDiff[ipoint+ipoint]-meanT)>rmsT*outlierCut){
123 : continue;
124 : }
125 0 : Double_t w= TMath::Gaus(dpoint);
126 0 : sumM+=magFull[ipoint+dpoint]*w;
127 0 : sumW+=w;
128 0 : }
129 0 : sumM/=sumW;
130 0 : reMod[ipoint]=sumM*TMath::Cos(phaseFull[ipoint]);
131 0 : imMod[ipoint]=sumM*TMath::Sin(phaseFull[ipoint]);
132 0 : }
133 : }
134 : //
135 : // 3.) Make a back FFT transformation
136 : //
137 0 : TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &fftLength, "C2R M K");
138 0 : fft_back->SetPointsComplex(reMod.GetMatrixArray(),imMod.GetMatrixArray());
139 0 : fft_back->Transform();
140 : TH1D *hb = 0;
141 0 : hb = (TH1D*)(TH1::TransformHisto(fft_back,hb,"Re"));
142 0 : hb->SetTitle("The backward transform result");
143 0 : hb->Scale(1./Double_t(fftLength));
144 0 : TGraph *grOutput = new TGraph(grInput);
145 0 : for (Int_t i=0; i<fftLength; i++){
146 : //hinput->SetBinContent(i, hb->GetBinContent(i-firstBin));
147 0 : grOutput->GetY()[i]=hb->GetBinContent(i+1);
148 : }
149 : //
150 : // 4.) in case specified streamer dump intermediat results for the checking purposes
151 : //
152 0 : if (pcstream!=NULL){
153 0 : (*pcstream)<<"fft"<<
154 0 : "hinput."<<hinput<< // input histogram
155 0 : "grInput.="<<&grInput<< // input graph
156 0 : "grOutput.="<<grOutput<< // output graph with removed outlier frequencies
157 : //
158 0 : "reFull.="<<&reFull<< // fft real part for original graph
159 0 : "imFull.="<<&imFull<< // fft imaginary part for original graph
160 0 : "reMod.="<<&reMod<<
161 0 : "imMod.="<<&imMod<<
162 : //
163 0 : "magFull.="<<&magFull<<
164 0 : "magMed.="<<&magMed<<
165 0 : "vecDiff.="<<&vecDiff<<
166 : //
167 0 : "meanT="<<meanT<< // mean difference
168 0 : "rmsT="<<rmsT<< // rms difference
169 : "\n";
170 : }
171 0 : delete hb;
172 : return grOutput;
173 0 : }
174 :
175 :
|
