TSpectrum and TSpectrum2

Dear ROOT wizards,

I am having problems with TSpectrum and TSpectrum2. About 1% of the time TSpectrum is not finding visible peaks. Also, TSpectrum2 doesn’t ever seem to work. I am using it in compiled code. I attach a snippet of the call:

TSpectrum* s = new TSpectrum();
s->Search(_histoRadius,2.," ",0.5);
TSpectrum2* s2 = new TSpectrum2();
s2->Search(_histoRadiusDCA);

please look in the checkhits directory and take a look at event 39. I see an obvious peak. Also, I have run 1000 events, found 9 bad ones total, and don’t see a clear pattern to tell me what characterizes these events.

In the TSpectrum2 calls, look at (in checkhits) at any of the RadiusDCA 2D histos: e.g., _histoRadiusDCA00002. It finds lots of peaks, none of them correct. The histos are created (using tfs) like this:

TH2D* _histoRadiusDCA = tfs->make<TH2D>(radiusDCAHisto.name(),radiusDCAHisto.title(),
					80,100.,500.,100,0.,100.);

I’ve tried TH2F and get the same results. One can see obvious peaks in the histogram. I did all this on one machine (ilcsim2.fnal.gov) and then decided to check another machine, so I also ran peak2.C on the original machine and fnalu (flxi05.fnal.gov for specificity) and the results look somewhat sensible, but there seem to be lots of ghosts in the upper right corner of the plots in all cases on both machines. I would think the ghosts should be randomly distributed.

Can someone help me out? Perhaps I need to tweak an argument or I’m doing something obviously wrong. Phillipe Canal is a couple of floors down so that might be a way to get quick info out of me.

–thanks, Robert Bernstein
houghtest_03.root (198 KB)

I do not see any problem with the 1D case if you set sigma=1 and a threshold >0.2
The same settings will help too for the 2D case. I am currently busy with many other things, but will look into TSpectrum2 in case of your 2D histos.

Rene

Dear Rene,

Thanks. It works interactively when I read in the histo but not when compiled for either set of choices in the search. I have seen similar problems with fits using MINUIT that don’t work when ROOT is compiled, but work perfectly when one writes out the histogram and runs ROOT interactively.

However, the 2D case fails when run either compiled or interactive. For the 1D case, I can live with the 0.9% failure rate and can probably tweak things. However, this is a wonderful add-on and when you can get to the 2D case I would appreciate it.

–thanks, Robert Bernstein

Hi Robert,

It would be a good help if you could provide a concrete piece of code showing the different behavior between compiled and interpreted code.
I mean for TSpectrum and also TMinuit.

Rene

Deari Rene,

I’m not entirely sure what you mean. The histogram I sent in the first e-mail (in the attachment) had no polymarker over the peak, so TSpectrum did not find it. The rest of the histos do have the polymarker. I sent the code snippet – what do you need me to send you?

On MINUIT, I realize I haven’t given you enough information to do anything useful. After I understand what you need, I can recreate the problem for TMinuit and document it adequately.

–thanks, Robert

Example of code that works in interpreted or compiled mode

Rene

#include "TFile.h" #include "TH1.h" #include "TSpectrum.h" void pic(Int_t event) { TFile *f= TFile::Open("houghtest_03.root"); if (event<10) h = (TH1F*)f->Get(Form("checkhits/_histoRadius0000%d",event)); else h = (TH1F*)f->Get(Form("checkhits/_histoRadius000%d",event)); TSpectrum s(50); h->Draw(); s.Search(h,1.," ",0.5); }

Dear Rene,

I noted you supplied a value of maxPositions (TSpectrum(50) )which I left as a default, i.e. TSpectrum(). When I set it to 50 (or any big number) it immediately worked and found the peak. I looked in the source and if maxPositions <= 0, the code sets maxPositions to 1. Now it does find multiple peaks in other events of mine even if you let it be the default at the constructor, so it’s seems that it’s some subtlety of the algorithm where if you start with too small a value of maxPositions, sometimes it finds none.

How did you know to do that? Should the code require the user to set a value?

I’m still baffled by the TSpectrum2 class but I will try explicit calls. However, the results of peaks2.C still looks odd to me.

–thanks so much, Robert