AliRoot Core  837871f (837871f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFFTsmoother.cxx
Go to the documentation of this file.
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  **************************************************************************/
16 
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 {
56  TGraph grInput(hinput);
57  Int_t fftLength = hinput->GetNbinsX();
58  TH1D *hm = 0;
59  TVirtualFFT::SetTransform(0);
60  hm = (TH1D*)hinput->FFT(hm, "MAG");
61  hm->SetTitle("Magnitude of the 1st transform");
62 
63  //
64  // 1.) Make FFT of the input histogram
65  //
66  TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
67  TVectorD reFull(fftLength);
68  TVectorD imFull(fftLength);
69  TVectorD magFull(fftLength);
70  TVectorD magS(fftLength);
71  TVectorD phaseFull(fftLength);
72  fft->GetPointsComplex(reFull.GetMatrixArray(),imFull.GetMatrixArray());
73  for (Int_t ipoint=1; ipoint<fftLength/2-1; ipoint++){
74  magFull[ipoint]=TMath::Sqrt(reFull[ipoint]*reFull[ipoint]+imFull[ipoint]*imFull[ipoint]);
75  phaseFull[ipoint]=TMath::ATan2(imFull[ipoint], reFull[ipoint]);
76  }
77  //
78  TVectorD reMod(fftLength);
79  TVectorD imMod(fftLength);
80  TVectorD vecDiff(fftLength);
81  //
82  // 2.) identify and replace outlier frequencies
83  //
84  TVectorD magMed(fftLength);
85  for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
86  Int_t dedge=TMath::Min(ipoint, fftLength/2-ipoint-1);
87  Int_t window= TMath::Min(dedge, medianRange);
88  magMed[ipoint]=TMath::Median(1+2*window, &(magFull.GetMatrixArray()[ipoint-window]));
89  vecDiff[ipoint]=0;
90  if (magMed[ipoint]>0){
91  vecDiff[ipoint]=(magFull[ipoint]-magMed[ipoint]);
92  }
93  }
94  for (Int_t ipoint=1; ipoint<fftLength/2-1; ipoint++){
95  vecDiff[ipoint]= (magFull[ipoint]-(magMed[ipoint-1]+magMed[ipoint+1])*0.5);
96  }
97  //
98  Double_t meanT, rmsT;
99  TStatToolkit::EvaluateUni(fftLength/2,vecDiff.GetMatrixArray(), meanT, rmsT, 0.95*(fftLength/2));
100  if (smoothSigma<=0){
101  for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
102  reMod[ipoint]=reFull[ipoint];
103  imMod[ipoint]=imFull[ipoint];
104  if (ipoint<lowBand) continue;
105  if (ipoint>=fftLength/2-lowBand) continue;
106  if (TMath::Abs(vecDiff[ipoint]-meanT)>rmsT*outlierCut){
107  reMod[ipoint]=magMed[ipoint]*TMath::Cos(phaseFull[ipoint]);
108  imMod[ipoint]=magMed[ipoint]*TMath::Sin(phaseFull[ipoint]);
109  }
110  }
111  }else{
112  for (Int_t ipoint=0; ipoint<fftLength/2-1; ipoint++){
113  reMod[ipoint]=reFull[ipoint];
114  imMod[ipoint]=imFull[ipoint];
115  if (ipoint<lowBand) continue;
116  if (ipoint>=fftLength/2-lowBand) continue;
117  //
118  Double_t sumM=0, sumW=0;
119  for (Int_t dpoint=-4.*smoothSigma; dpoint<4.*smoothSigma; dpoint++){
120  if (dpoint+ipoint<0) continue;
121  if (dpoint+ipoint>fftLength/2-1) continue;
122  if (TMath::Abs(vecDiff[ipoint+ipoint]-meanT)>rmsT*outlierCut){
123  continue;
124  }
125  Double_t w= TMath::Gaus(dpoint);
126  sumM+=magFull[ipoint+dpoint]*w;
127  sumW+=w;
128  }
129  sumM/=sumW;
130  reMod[ipoint]=sumM*TMath::Cos(phaseFull[ipoint]);
131  imMod[ipoint]=sumM*TMath::Sin(phaseFull[ipoint]);
132  }
133  }
134  //
135  // 3.) Make a back FFT transformation
136  //
137  TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &fftLength, "C2R M K");
138  fft_back->SetPointsComplex(reMod.GetMatrixArray(),imMod.GetMatrixArray());
139  fft_back->Transform();
140  TH1D *hb = 0;
141  hb = (TH1D*)(TH1::TransformHisto(fft_back,hb,"Re"));
142  hb->SetTitle("The backward transform result");
143  hb->Scale(1./Double_t(fftLength));
144  TGraph *grOutput = new TGraph(grInput);
145  for (Int_t i=0; i<fftLength; i++){
146  //hinput->SetBinContent(i, hb->GetBinContent(i-firstBin));
147  grOutput->GetY()[i]=hb->GetBinContent(i+1);
148  }
149  //
150  // 4.) in case specified streamer dump intermediat results for the checking purposes
151  //
152  if (pcstream!=NULL){
153  (*pcstream)<<"fft"<<
154  "hinput."<<hinput<< // input histogram
155  "grInput.="<<&grInput<< // input graph
156  "grOutput.="<<grOutput<< // output graph with removed outlier frequencies
157  //
158  "reFull.="<<&reFull<< // fft real part for original graph
159  "imFull.="<<&imFull<< // fft imaginary part for original graph
160  "reMod.="<<&reMod<<
161  "imMod.="<<&imMod<<
162  //
163  "magFull.="<<&magFull<<
164  "magMed.="<<&magMed<<
165  "vecDiff.="<<&vecDiff<<
166  //
167  "meanT="<<meanT<< // mean difference
168  "rmsT="<<rmsT<< // rms difference
169  "\n";
170  }
171  delete hb;
172  return grOutput;
173 }
174 
175 
TTreeSRedirector * pcstream
void EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
TGraph * ReplaceOutlierFrequenciesMedian(TH1 *hinput, Double_t outlierCut, Int_t medianRange, Float_t smoothSigma, Int_t lowBand, TTreeSRedirector *pcstream)
class TVectorT< Double_t > TVectorD