Wrapper function with TSpectrum

NgJunKai · April 19, 2023, 8:19am

i wrote a macro as in below

void identify_hist_peaks2(TH1D* hist,int maxpositions=3, double resolution=1.0, bool draw=false, bool addToHist=false,int decimal_places=2,string unit="MHz"){

	cout << unit << endl;

    TSpectrum * peak_finder = new TSpectrum(maxpositions, resolution); // (int maxpositions, double resolution = 1)
    int npeaks = peak_finder->Search(hist,2.0,"",0.1); // int TSpectrum::Search(const TH1* hist, double sigma = 2, const char* option = "", double threshold = 0.050000000000000003)
    cout << "found " << npeaks << " peaks." << endl;
    TPolyMarker* pm = (TPolyMarker*) hist->FindObject("TPolyMarker");
    
    if(addToHist){hist->GetListOfFunctions()->Add(pm->Clone()); };
    if(draw){pm->Draw(); };
     // caveat: if not draw, the TPoly will resides in the ListOfFunction of that histo.
     //         if draw, the TPoly need to be attached to the ListOfFunction, in order to be drawn

    for (int pk=0; pk<npeaks;pk++){
        TLatex* text = new TLatex(pm->GetX()[pk]*1.05, pm->GetY()[pk]*0.9, Form("%.*f %s", decimal_places, pm->GetX()[pk],unit.c_str()));
        text->SetTextColor(2);
        text->SetTextSize(0.04);
        if(draw){text->DrawClone(); };
        if(addToHist){hist->GetListOfFunctions()->Add(text); };
    };
}

then in pyroot, i run ROOT.gInterpreter.Declare('# include "test.h"') (also tried ROOT.gROOT.ProcessLine()). I have some fft histograms that i use the ‘imported’ function to find and label the peaks, below is the code that does the job.

c = r.TCanvas("cc","c",800,400)
hist_fft.Draw()
ROOT.identify_hist_peaks2(hist_fft,maxpositions=3, resolution = 1.0,draw=True,unit="zzz") 
c.Draw()

and wierd things happens:
if i do exactly the same as in the above code, i got error saying could not convert argument 4 (tuple assignment index out of range).
on the other hand, if i keep only the pass only the hist and draw=True, then i got (jupyter) kernel restart.
is my code buggy or i just missed out something ?

ROOT Version: 6.24 (conda-forge)
Platform: ubuntu 20.04(wsl)
Compiler: gcc9

couet · April 19, 2023, 8:49am

As your macro is a C++ script, have you tried to use it first in plain C++ interpreter ?

NgJunKai · April 19, 2023, 8:52am

not yet. will try …

NgJunKai · April 19, 2023, 9:26am

i prepare 3 files

the problematic file test.h (full content)

void identify_hist_peaks2(TH1D* hist,int maxpositions=3, double resolution=1.0, bool draw=false, bool addToHist=false,int decimal_places=2,string unit="MHz"){

	cout << unit << endl;

    TSpectrum * peak_finder = new TSpectrum(maxpositions, resolution); // (int maxpositions, double resolution = 1)
    int npeaks = peak_finder->Search(hist,2.0,"",0.1); // int TSpectrum::Search(const TH1* hist, double sigma = 2, const char* option = "", double threshold = 0.050000000000000003)
    cout << "found " << npeaks << " peaks." << endl;
    TPolyMarker* pm = (TPolyMarker*) hist->FindObject("TPolyMarker");
    
    if(addToHist){hist->GetListOfFunctions()->Add(pm->Clone()); };
    if(draw){pm->Draw(); };
     // caveat: if not draw, the TPoly will resides in the ListOfFunction of that histo.
     //         if draw, the TPoly need to be attached to the ListOfFunction, in order to be drawn

    // int npeaks = pm->GetN();
    for (int pk=0; pk<npeaks;pk++){
        //cout << pm->GetX()[pk] << endl;
        //TString formatted_text = Form("%.*f %s", decimal_places, pm->GetX()[pk],unit.c_str());
        TLatex* text = new TLatex(pm->GetX()[pk]*1.05, pm->GetY()[pk]*0.9, Form("%.*f %s", decimal_places, pm->GetX()[pk],unit.c_str()));
        text->SetTextColor(2);
        text->SetTextSize(0.04);
        if(draw){text->DrawClone(); };
        if(addToHist){hist->GetListOfFunctions()->Add(text); };
    };
}

the file used to do fft ToolsCollection.h

#ifndef TOOLSCOLLECTION_H
#define TOOLSCOLLECTION_H

#include <iostream>

using namespace std;

TH1D* fft(TH1D* hist, string fftname="default_fft_name", string given_title="FFT", double fstart=0.01){
    double bin_width = hist->GetBinWidth(1);
    cout << bin_width << endl;
    
    TString title = Form("%s ;f [MHz]; FFT Magnitude [arb.]",given_title.c_str());
    
    TH1D* fft_hist = new TH1D(fftname.c_str(), title.Data(), hist->GetNbinsX(), 0, 1.0/bin_width);
    hist->FFT(fft_hist, "MAG");
    fft_hist->GetXaxis()->SetRangeUser(fstart,fft_hist->GetXaxis()->GetXmax()/2); // fstart to get rid of the DC component
    return fft_hist;
}



#endif

file that takes all code together testtest.C

#include <iostream>
#include <string>

#include "root_tools/ToolsCollection.h"
#include "test.h"

using namespace std;

void testtest(){
    
int npts = 4000;

TF1* func = new TF1("func","1 + 0.5*cos(14.48*x + 0)",0,npts*0.1492); 
func->SetNpx(npts*2);

TH1D* hist = new TH1D("hist","hist",npts,0,npts*0.1492);
    
int N = 100000;
for(int i=0;i<N;i++){
    double t = func->GetRandom();
    hist->Fill(t);
    };   

TCanvas c = TCanvas("cc","c",400,400);
TH1D* hist_fft = fft(hist,"fft_name","FFT",0.01); // put here xwork !!!
hist_fft->Draw();
identify_hist_peaks2(hist_fft,true,"GHz");
c.Draw();
c.Print("test.png");
exit(0);
}

let say i only want to change the two argument, ‘draw=trueandunit="GHz’. when i did this, i got error: no matching function for call to 'identify_hist_peaks2'

yus · April 19, 2023, 9:53am

Hi,

this can’t work. The function has 7 arguments. By providing only three you try to change the first three, but, for example, the second one is supposed to be an int, while you pass true to it. You should do the following instead:

identify_hist_peaks2(hist_fft,3, 1.0, true, false, 2, "GHz");

NgJunKai · April 19, 2023, 10:59am

thanks for pointing that out.
so if i want to change the unit, it is better to do in this way.

hist_fft->Draw();
int maxpositions=3;
double resolution=1.0;
bool draw=true;
bool addToHist=false;
int decimal_places=2;
string unit="GHz";
identify_hist_peaks2(hist_fft,maxpositions,resolution,draw,addToHist,decimal_places,unit);

or if i change this frequently, it is better to put unit in the from at the time of definition, so i can safely change only the first few without touching the rest.

sorry i code frequently in python (and less in c++). so i just mixed up both. is it almost the same rules (see image below) applied to c++ ?

yus · April 19, 2023, 11:34am

NgJunKai:

so if i want to change the unit, it is better to do in this way.

hist_fft->Draw();
int maxpositions=3;
double resolution=1.0;
bool draw=true;
bool addToHist=false;
int decimal_places=2;
string unit="GHz";
identify_hist_peaks2(hist_fft,maxpositions,resolution,draw,addToHist,decimal_places,unit);

yes, or simply

identify_hist_peaks2(hist_fft, 3, 1., true, false, 2, "GHz");

In C++, the order is very important. You should either provide all arguments (in the same order these are used in the function) or only up to the last one you want to change from the default value (then the rest of the arguments will use the default values if provided; if not, it’s a compilation error; the order of the first N passed values should match the order of the corresponding arguments of the function).

NgJunKai · April 19, 2023, 11:55am

ok got it. to put it simply, there is no “keyword argument” in c++

yus · April 19, 2023, 12:01pm

Exactly. Though there are several tricks to introduce keyword arguments in C++ (see e.g. Named Arguments in C++ - Fluent C++), there is no default support for these.

system · May 3, 2023, 12:01pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.