Problem with peak search tool of TSpectrum

ROOT
1

Good afternoon,

I am analyzing a spectrum using the tools provided by the TSpectrum class.
There is some bug in my macro I could not find yet.

When performing a background subtraction from the original histogram

 // Calculating and substracting background
      //
  TSpectrum* spect= new TSpectrum(12,1.) ;
  TSpectrum* spectsb= new TSpectrum(30,2.0) ;
  hdyne1back=(TH1F*)spect->Background(hdyne1,160,"BackDecreasingWindow BackOrder8 nosmoothing Compton same") ;
  hdyne1back->SetLineColor(kGreen) ; 
  hdyne1sb->Add(hdyne1,hdyne1back,1,-1) ;
  hdyne1sb->SetLineColor(kBlue) ;

, the result looks fine.

Then I want to search peaks on the histogram obtained after subtracting the background:

// Peak Research
       
  spectsb->Search(hdyne1sb,2.5,"nobackground",0.2) ;
  TList *functions = hdyne1sb->GetListOfFunctions();
  TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
  pm->SetMarkerSize(0.7) ;
  pm->SetMarkerColor(kRed) ;
  npeakssb=spectsb->GetNPeaks() ;
  cout << " NPeaks= " << npeakssb << endl ;
  xpeaksb=spectsb->GetPositionX() ;
  xpeakmaxsb=0.0 ; // Peak Max (for 137Cs)

  //Reordering peaks list
  for(Int_t ipc=0;ipc<npeakssb;ipc++) {
    cout << "Dyne1 X Peak " << ipc << " " << xpeaksb[ipc] << endl ;
    if(xpeaksb[ipc]>=xpeakmaxsb) xpeakmaxsb=(Double_t)xpeaksb[ipc] ;
    for(Int_t ipcorder=ipc+1;ipcorder<npeakssb;ipcorder++) { // Loop of reordering
	  

      if(xpeaksb[ipcorder]<xpeaksb[ipc]) { // Begin if
	a=xpeaksb[ipc] ;
	b=xpeaksb[ipcorder] ;

	xpeaksb[ipcorder]=a ;
	xpeaksb[ipc]=b ;

      } // Endif
      
    } // End loop of reordering 
  }
  cout << "Max peak: " << xpeakmaxsb << endl ;

The fact is that the Polymarker positions along the X-axis are not consistent with the values stored in the xpeaksb pointer. Indeed, the peak located around channel 600 is not stored although a polymarker is drawn on the histogram at this position … When dumping on screen I have:

Reading File: …/data09122011/22Na_60Co_137Cs_LaBr3_2x2_R7723-100_1800V.dat.root
Entries: 100001
NPeaks= 4
Dyne1 X Peak 0 771.123
Dyne1 X Peak 1 771.123
Dyne1 X Peak 2 1353.84
Dyne1 X Peak 3 1528.75
Max peak: 1528.75

The second peak is repeted twice …

Up to now I did not find out where the bug could be. Does anyone have any suggestion ? I am using Root version 5.27/06b.

I attached the root file and the macro.
In the macro you could have to change the path to the root file to execute the macro properly:

 char* filename="../data09122011/22Na_60Co_137Cs_LaBr3_2x2_R7723-100_1800V.dat.root" ;

Thanks a lot for the help

2

Just to make sure we have the same result, here is the picture I get with your macro. Do you get the same one ?

ca3.gif696×572 18.3 KB

3

Yes that’s correct. I am sorry I just found out that this post was submitted in the wrong forum.

I made another investigation by reading the X positions of polymarkers folowing:

spectsb->Search(hdyne1sb,3.0,"nobackground",0.2) ;
  TList *functions = hdyne1sb->GetListOfFunctions();
  TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
  pm->SetMarkerSize(0.7) ;
  pm->SetMarkerColor(kRed) ;
  cout << "Number of polymarkers :" << pm->GetN() << endl ;
  Double_t* pmarrayX=pm->GetX() ; // Find polymarkers position
  npeakssb=spectsb->GetNPeaks() ;

and then to print on screen:

//Reordering peaks list
  for(Int_t ipc=0;ipc<npeakssb;ipc++) {
    cout << "Dyne1 X Peak " << ipc << " " << xpeaksb[ipc] << endl ;
   cout << "Dyne1 X Polymarker " << ipc << " " << pmarrayX[ipc] << endl ; // print Polymarkers position
    if(xpeaksb[ipc]>=xpeakmaxsb) xpeakmaxsb=(Double_t)xpeaksb[ipc] ; 
    for(Int_t ipcorder=ipc+1;ipcorder<npeakssb;ipcorder++) { // Loop of reordering
	  

      if(xpeaksb[ipcorder]<xpeaksb[ipc]) { // Begin if
	a=xpeaksb[ipc] ;
	b=xpeaksb[ipcorder] ;

	xpeaksb[ipcorder]=a ;
	xpeaksb[ipc]=b ;

      } // Endif
      
    } // End loop of reordering 
  }
  cout << "Max peak: " << xpeakmaxsb << endl ;

The result is the following on screen:

Reading File: …/data09122011/22Na_60Co_137Cs_LaBr3_2x2_R7723-100_1800V.dat.root
Entries: 100001
Number of polymarkers :4
NPeaks= 4
Dyne1 X Peak 0 771.123
Dyne1 X Polymarker 0 771.123
Dyne1 X Peak 1 771.123
Dyne1 X Polymarker 1 599.207
Dyne1 X Peak 2 1353.84
Dyne1 X Polymarker 2 1353.84
Dyne1 X Peak 3 1530.75
Dyne1 X Polymarker 3 1530.75
Max peak: 1530.75

When using the polymarkers coordinates I found the correct peak positions !

The result is the same when using Root version 5.32/00 .

4 
Dyne1 X Peak 0 771.123
Dyne1 X Polymarker 0 771.123
Dyne1 X Peak 1 771.123
Dyne1 X Polymarker 1 599.207
Dyne1 X Peak 2 1353.84
Dyne1 X Polymarker 2 1353.84
Dyne1 X Peak 3 1530.75
Dyne1 X Polymarker 3 1530.75

So according to this output the peak postion and the marker position are the same … so that’s correct, right ? … where is the problem in that case ?

5

Now moved in the right one “ROOT support”

6

The problem is why polymarkers coordinates for the peaks are different from the values stored in the xpeaks pointer … ?

7

I am lost … they look the same on the output given by your macro…
I think I have not understood your point.

8

Hello,

I printed the data stored in the polymarker X position pointer (Dyne1 X Polymarker) and the pointer that stores data from the Search function (Dyne1 X Peak). There are not the same because in the second case, the peak located at channel 600 is not found, whereas the peak at channel 770 is quoted twice

Dyne1 X Peak 0 771.123
Dyne1 X Peak 1 771.123
Dyne1 X Peak 2 1353.84
Dyne1 X Peak 3 1530.75

Dyne1 X Polymarker 0 771.123
Dyne1 X Polymarker 1 599.20
Dyne1 X Polymarker 2 1353.84
Dyne1 X Polymarker 3 1530.75

9

Oh I see now… and which one do you think is right ? (I try to understand because I am not the author of that code).

10

The results stored in the polymarker look fine (you can see small red markers on top of the peaks). What I don’t understand is why in the first list the first peak does not appear and the second is repeted twice …
I expected that results are the same for both pointers

Could it be some memory side effect in the implementation in the Search function of the TSpectrum class ?

Even if I can proceed further, I would like to understand such difference, to be sure there is no hidden bug …

11

I think your code is faulty.

I get:

If I simplify the the loop to reorder peaks as follow:

  //Reordering peaks list
  for(Int_t ipc=0;ipc<npeakssb;ipc++) {
    cout << "Dyne1 X Peak " << ipc << " " << xpeaksb[ipc] << endl ;
     /*
    if(xpeaksb[ipc]>=xpeakmaxsb) xpeakmaxsb=(Double_t)xpeaksb[ipc] ;
    for(Int_t ipcorder=ipc+1;ipcorder<npeakssb;ipcorder++) { // Loop of reordering
	  

      if(xpeaksb[ipcorder]<xpeaksb[ipc]) { // Begin if
	a=xpeaksb[ipc] ;
	b=xpeaksb[ipcorder] ;

	xpeaksb[ipcorder]=a ;
	xpeaksb[ipc]=b ;

      } // Endif
      
    } // End loop of reordering */
  }
12

You’re right ! Small error, big consequences. Thanks a lot !!