TestChannel.C
The following ROOT macro plots the input, output and reconstructed signals of a single GET channel:
- the input sample has 2 Gauss component
- some input noise is added to the signal (GETNoise class)
- the response function is analystical (GETResponseFctAget class)
- some noise is added to the output signal (GETNoise class)
- the reconstruction is performed with filtering (GETFilterAgetGaus class)
The TestChannel.C example requires the functions defined in the TestCore.C. In this example, it is necessary to compile the ROOT macro files, because they use function pointers, that are not supported in the ROOT command interpreter.
//======================================================================
// ROOT macros channel test
//
//
//
//======================================================================
/*
These macros require to load the GET libraries
gROOT->Macro ( (getenv("GET_BASE")+string("/root/GETRootLibs.C")).c_str() );
To load the macro file
.L TestCore.C+
.L TestChannel.C+
TestChannel();
*/
//======================================================================
#include "TestCore.H"
void TestChannel ( )
{
//------------------------------------------------------------
u_int n = 256; // sample size
double dt = 0.02; // sampling time (us)
double flin = 10.0L; // input fluctuation (uA)
double flout = 2.0L; // output fluctuation
double fcut = 2.5L; // filter frequency cut (MHz)
double twid = 0.25L; // signal width (FWHM - us)
double peak = 0.40L; // response peaking time (us)
double tmin = -0.5*dt;
double tmax = (n-0.5)*dt;
double df = 0.5L/(dt*n);
double fmin = - 0.5*df;
double fmax = 0.5L/dt + 0.5*df;
double resp_fact = 20.L; // graph scale factors
double filt_fact = 10.L;
double fref_fact = 2000.L;
Color_t col_in0 = kRed+1;
Color_t col_in = kOrange-3;
Color_t col_out0 = kGreen+3;
Color_t col_out1 = kGreen+2;
Color_t col_out = kGreen+1;
Color_t col_rec = kRed;
Color_t col_resp = kBlue+1;
Color_t col_filt = 1;
//------------------------------------------------------------
// Create GET channel and signals
cerr << endl;
cerr << "* Creating GET channel" << endl;
GETChannel channel(n,dt);
// response function
cerr << " - response function" << endl;
GETResponseFctAget resp ( n, dt, 50.L, peak, 0.L, 0.L );
channel.SetResponse ( resp );
// create pure input and output signals (not in channel: for
// display only
// Charge Center FWHM
// fC fC us us us
double SinPar[5] = { 10., 5., 1.5, 2.0, twid };
cerr << " - setting pure input signal" << endl;
GETSample input0 ( n, dt );
input0.SetFunctionData ( FctGaus2, SinPar );
cerr << " N = " << input0.GetDimension() << endl;
cerr << " dT = " << input0.GetTimeStep() << endl;
cerr << " T[0] = " << input0.GetTimeMin() << endl;
cerr << " T[N] = " << input0.GetTimeEnd() << endl;
cerr << " dF = " << input0.GetFrequencyStep() << endl;
cerr << " Fc = " << input0.GetCriticalFrequency() << endl;
// pure output signal
cerr << " - generating pure output signal" << endl;
GETSample output0 ( n, dt );
resp.ProcessSample ( input0, output0, channel );
// input noise generator
// it is not set to the channel processes, since it is
// applied externally, for comparison of results without
// or with input noise
cerr << " - input noise generator" << endl;
cerr << " fluctuation FWHM = " << flin << " uA" << endl;
GETNoise noise_in ( n, dt, flin );
////channel.SetNoiseIn ( noise_in );
// noisy input signal
GETSample input1 ( input0 );
cerr << " - generating channel input with noise" << endl;
noise_in.ProcessSample ( input0, input1, channel );
cerr << " FFT(In) = " << channel.TstSample().IsTransformSet() << endl;
channel.TstSample() = input1;
// output signal with input noise
cerr << " - generating output signal with input noise" << endl;
GETSample output1 ( n, dt );
resp.ProcessSample ( channel.TstSample(), output1, channel );
// output noise generator
cerr << " - output noise generator" << endl;
cerr << " fluctuation FWHM = " << flout << " uA" << endl;
GETNoise noise_out ( n, dt, flout );
channel.SetNoiseOut ( noise_out );
// processing filter
cerr << " - processing filter" << endl;
GETFilterGaus filt ( n, dt, fcut, twid );
//GETFilterAget filt ( n, dt, fcut, peak );
////GETFilterAgetGaus filt ( n, dt, fcut, twid, peak );
////GETFilterLowN filt1 ( n, dt, 4, fcut );
////GETFilterLowN filt2 ( n, dt, 4, -fcut );
////GETFilter filt ( filt1, filt2 );
channel.SetFilter ( filt );
// output signal
cerr << " - output signal" << endl;
channel.ProcessOutFromTst();
cerr << " fct(In) = " << channel.TstSample().IsFunctionSet() << endl;
cerr << " FFT(In) = " << channel.TstSample().IsTransformSet() << endl;
cerr << " FFT(Out) = " << channel.OutSample().IsTransformSet() << endl;
// reconstructed signal
cerr << " - reconstructed signal" << endl;
channel.ProcessRecFromOut();
//------------------------------------------------------------
// Create graphs for various signals
cerr << endl;
cerr << "* Creating graphs" << endl;
// response function
TGraph * SrespGr = resp.GetResponse().CreateFunctionGraph(resp_fact);
TGraph * PrespGr = resp.GetResponse().CreatePowerGraph(resp_fact);
SrespGr->SetLineColor ( col_resp );
SrespGr->SetLineStyle ( 1 );
SrespGr->SetLineWidth ( 2 );
PrespGr->SetLineColor ( col_resp );
PrespGr->SetLineStyle ( 1 );
PrespGr->SetLineWidth ( 2 );
// filter
TGraph * SfiltGr = filt.GetFilter().CreateFunctionGraph(filt_fact);
TGraph * PfiltGr = filt.GetFilter().CreatePowerGraph(filt_fact);
SfiltGr->SetLineColor ( col_filt );
SfiltGr->SetLineStyle ( 1 );
SfiltGr->SetLineWidth ( 2 );
PfiltGr->SetLineColor ( col_filt );
PfiltGr->SetLineStyle ( 1 );
PfiltGr->SetLineWidth ( 2 );
// filter reference function
/*
TGraph * SfrefGr = filt.GetReference().CreateFunctionGraph(fref_fact);
TGraph * PfrefGr = filt.GetReference().CreatePowerGraph(fref_fact);
TGraph * SfrefGr = filt.GetReference().CreateFunctionGraph(fref_fact);
TGraph * PfrefGr = filt.GetReference().CreatePowerGraph(fref_fact);
SfrefGr->SetLineColor ( col_filt );
SfrefGr->SetLineStyle ( 2 );
SfrefGr->SetLineWidth ( 2 );
PfrefGr->SetLineColor ( col_filt );
PfrefGr->SetLineStyle ( 2);
PfrefGr->SetLineWidth ( 2 );
*/
// pure input
TGraph * Sin0Gr = input0.CreateFunctionGraph();
TGraph * Pin0Gr = input0.CreatePowerGraph();
Sin0Gr->SetLineColor ( col_in0 );
Sin0Gr->SetLineStyle ( 2 );
Sin0Gr->SetLineWidth ( 3 );
Pin0Gr->SetLineColor ( col_in0 );
Pin0Gr->SetLineStyle ( 2 );
Pin0Gr->SetLineWidth ( 3 );
// noisy input
TGraph * SinGr = channel.TstSample().CreateFunctionGraph();
TGraph * PinGr = channel.TstSample().CreatePowerGraph();
SinGr->SetLineColor ( col_in );
SinGr->SetLineStyle ( 1 );
SinGr->SetLineWidth ( 2 );
PinGr->SetLineColor ( col_in );
PinGr->SetLineStyle ( 1 );
PinGr->SetLineWidth ( 2 );
// pure output
TGraph * Sout0Gr = output0.CreateFunctionGraph();
TGraph * Pout0Gr = output0.CreatePowerGraph();
Sout0Gr->SetLineColor ( col_out0 );
Sout0Gr->SetLineStyle ( 2 );
Sout0Gr->SetLineWidth ( 3 );
Pout0Gr->SetLineColor ( col_out0 );
Pout0Gr->SetLineStyle ( 2 );
Pout0Gr->SetLineWidth ( 3 );
// output with input noise
TGraph * Sout1Gr = output1.CreateFunctionGraph();
TGraph * Pout1Gr = output1.CreatePowerGraph();
Sout1Gr->SetLineColor ( col_out1 );
Sout1Gr->SetLineStyle ( 1 );
Sout1Gr->SetLineWidth ( 2 );
Pout1Gr->SetLineColor ( col_out1 );
Pout1Gr->SetLineStyle ( 1 );
Pout1Gr->SetLineWidth ( 2 );
// output signal
TGraph * SoutGr = output2.CreateFunctionGraph();
TGraph * PoutGr = output2.CreatePowerGraph();
////TGraph * SoutGr = channel.OutSample().CreateFunctionGraph();
////TGraph * PoutGr = channel.OutSample().CreatePowerGraph();
SoutGr->SetLineColor ( col_out );
SoutGr->SetLineStyle ( 1 );
SoutGr->SetLineWidth ( 2 );
PoutGr->SetLineColor ( col_out );
PoutGr->SetLineStyle ( 1 );
PoutGr->SetLineWidth ( 2 );
// reconstructed signal
TGraph * SrecGr = channel.RecSample().CreateFunctionGraph();
TGraph * PrecGr = channel.RecSample().CreatePowerGraph();
SrecGr->SetLineColor ( col_rec );
SrecGr->SetLineStyle ( 1 );
SrecGr->SetLineWidth ( 2 );
PrecGr->SetLineColor ( col_rec );
PrecGr->SetLineStyle ( 1 );
PrecGr->SetLineWidth ( 2 );
//------------------------------------------------------------
// Create plot
cerr << endl;
cerr << "* Creating Plots" << endl;
gStyle->SetPadLeftMargin (0.07);
gStyle->SetPadRightMargin (0.03);
gStyle->SetPadTopMargin (0.04);
gStyle->SetPadBottomMargin(0.12);
TCanvas * canvas = new TCanvas ( "TestChannel", "TestChannel", 1200, 900 );
canvas->Divide ( 1, 2 );
// - Signals plot
canvas->cd ( 1 );
TH1D * hsig = new TH1D ( "SignalFrame", "", n, tmin, tmax );
hsig->SetMinimum ( -9. );
hsig->SetMaximum ( 59. );
hsig->SetLineColor ( 0 );
hsig->SetLineWidth ( 0 );
hsig->GetXaxis()->SetTitle ( "time (#mus)" );
hsig->GetYaxis()->SetTitle ( "amplitude (fC/#mus = nA)" );
hsig->GetXaxis()->CenterTitle ( );
hsig->GetYaxis()->CenterTitle ( );
hsig->GetYaxis()->SetTitleOffset ( 0.6 );
hsig->GetXaxis()->SetTitleSize ( 0.05 );
hsig->GetYaxis()->SetTitleSize ( 0.05 );
hsig->GetYaxis()->SetTickLength ( 0.015 );
hsig->DrawCopy();
TLine * line = new TLine ( tmin, 0., tmax, 0. );
line->SetLineWidth(1);
line->SetLineColor(1);
line->SetLineStyle(1);
line->Draw();
SrespGr->Draw("L");
SfiltGr->Draw("L");
Sin0Gr->Draw("L");
SinGr->Draw("L");
Sout0Gr->Draw("L");
Sout1Gr->Draw("L");
SoutGr->Draw("L");
SrecGr->Draw("L");
// - PSD plot
canvas->cd ( 2 )->SetLogy();
TH1D * hpsd = new TH1D ( "PowerFrame", "", n, fmin, fmax );
hpsd->SetMinimum ( 5.e-10 );
hpsd->SetMaximum ( 2.e4 );
hpsd->SetLineColor ( 0 );
hpsd->SetLineWidth ( 0 );
hpsd->GetXaxis()->SetTitle ( "frequency (MHz)" );
hpsd->GetYaxis()->SetTitle ( "power spectral density" );
hpsd->GetXaxis()->CenterTitle ( );
hpsd->GetYaxis()->CenterTitle ( );
hpsd->GetYaxis()->SetTitleOffset ( 0.6 );
hpsd->GetXaxis()->SetTitleSize ( 0.05 );
hpsd->GetYaxis()->SetTitleSize ( 0.05 );
hpsd->GetYaxis()->SetNdivisions ( 7 );
hpsd->GetYaxis()->SetTickLength ( 0.015 );
hpsd->DrawCopy();
PrespGr->Draw("L");
PfiltGr->Draw("L");
//PfrefGr->Draw("L");
Pin0Gr->Draw("L");
PinGr->Draw("L");
Pout0Gr->Draw("L");
Pout1Gr->Draw("L");
PoutGr->Draw("L");
PrecGr->Draw("L");
//------------------------------------------------------------
// Legend
canvas->cd(1);
double leg_x = 0.68;
double leg_dy = 0.06;
TLegend * legend = new TLegend ( leg_x, 0.92-8*leg_dy, 0.95, 0.92 );
legend->SetFillColor(0);
legend->SetFillStyle(0);
legend->SetShadowColor(0);
legend->SetBorderSize(0);
legend->SetLineWidth(2);
legend->AddEntry ( PrespGr, ("AGet response function (x" + GGetString(resp_fact) + ")").c_str(), "l" );
legend->AddEntry ( PfiltGr, ("Filter (x" + GGetString(filt_fact) + ")").c_str(), "l" );
//legend->AddEntry ( PfrefGr, ("Filter ref. function (x" + GGetString(fref_fact) + ")").c_str(), "l" );
legend->AddEntry ( Pin0Gr, "Input signal (clean)", "l" );
legend->AddEntry ( PinGr, "Input signal (noisy)", "l" );
legend->AddEntry ( Pout0Gr, "Output signal (clean)", "l" );
legend->AddEntry ( Pout1Gr, "Output signal (input noise)", "l" );
legend->AddEntry ( PoutGr, "Output signal (in/out noise)", "l" );
legend->AddEntry ( PrecGr, "Reconstructed signal", "l" );
legend->Draw();
//------------------------------------------------------------
canvas->cd(1)->SaveAs("plots/TestChannelSignal.eps");
canvas->cd(2)->SaveAs("plots/TestChannelPower.eps");
//------------------------------------------------------------
// Cleaning
delete hsig;
delete hpsd;
}
//======================================================================
/*
*/
