Problematic Rebinning resulting in a spike in the middle

ROOT
1

I am trying to do something which is rather simple.
I am trying to rebin a histogram by a factor of 4.
However there is always a spike at a specific channel(512 which is the middle channel of the specified x axis range), like the one shown in the figure

1280×800

Just for reference the input histogram is the following

1280×800

My code that does the rebinning is the following

[code]#include “Riostream.h”
void rebin(char * file_c) {

//*****************************************************
//* First executes the script “evnt2dat”.
//* The script accepts as an input a .evnt file
//* and creates a .dat file.
//* Then root creates the ntuple, histogram from each detector
//* and the DE-E scatter plot
//*
//* Execute it using
//* root -l ‘ntuple.C(“filename”)’
//* Note that the extension .evnt MUST NOT be used
//* and evnt2dat script must be on the same directory as
//* this macro.
//*****************************************************

TString file(file_c);
//gSystem->Exec(TString::Format("./evnt2dat %s",file.Data()));//Executes the script evnt2dat
TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
dir.ReplaceAll(“rebin.C”,"");
dir.ReplaceAll("/./","/");
ifstream in;
in.open(TString::Format("%s%s.dat",dir.Data(),file.Data()));

gROOT->SetStyle(“Plain”);

Float_t x,y;
Int_t events = 0;
Int_t channels=4096;//Change it according to the acquisition channels
//TFile *f = new TFile(TString::Format("%s.root",file.Data()),“RECREATE”);
//Create the histograms
TH1F *histo = new TH1F(“hist”,“histo”,channels,1,channels);

//Create the ntuple
//TNtuple *ntuple = new TNtuple(“ntuple”,“DE/E analysis”,“x:y”);

//Fill the histos and the ntuple
while (1) {
in >> x >> y;
if (!in.good()) break;
histo->Fill(int(x),y);
//ntuple->Fill(de1,e1,de2,e2,de3,e3);
events++;
}
printf(“Found %d events\n”,events);

TCanvas *c1 = new TCanvas(“c1”, “c1”,225,219,700,530);
histo->Draw(“L”);
histo->SetBins(histo->GetNbinsX(), 0, 1024);
c1->Update();
histo->Draw();
histo->Rebin(4);
TString histfilename =TString::Format("%s_Rebinned.dat",file.Data());
SingleExportAscii(histo,histfilename);

gSystem->Exec(TString::Format("./copy %s_Rebinned.dat",file.Data()));
//gSystem->Exec(“gnuplot plot”);

in.close();
//f->Write();
}

/**

  • \brief Export Single Histogram into ASCII file
    /
    Bool_t SingleExportAscii(TH1     hist, TString &filename, TString folder="", TString separator="\t")
    {
    Int_t i,j;
    Double_t xcenter, xwidth;
    Bool_t success=kFALSE;
    //filename = folder + hist->GetName() + “.dat”;
    ofstream file_out(filename);
    file_out << "# Output " << hist->ClassName() << “: " << hist->GetName() << " (” << hist->GetTitle() << “)\n”;
    if (hist->GetDimension()==1)
    {
    file_out << “# BinCenter” << separator << “Content\n”;
    for (i=1; i<=hist->GetNbinsX(); i++)
    file_out << int(hist->GetBinCenter(i)-0.5) << separator << hist->GetBinContent(i) << endl;
    if (i>1)
    success=kTRUE;
    }
    else if (hist->GetDimension()==2)
    {
    file_out << “# xBinCenter” << separator << “yBinCenter” << separator << “Content” << separator << “xBinHalfWidth” << separator << “yBinHalfWidth” << separator << “Error” << endl;
    for (i=1; i <= hist->GetNbinsX(); i++)
    {
    xcenter = hist->GetXaxis()->GetBinCenter(i);
    xwidth = hist->GetXaxis()->GetBinWidth(i)/2;
    for (j=1; j <= hist->GetNbinsY(); j++)
    file_out << xcenter << separator << hist->GetYaxis()->GetBinCenter(j) << separator << hist->GetBinContent(i,j) << separator << xwidth << separator << hist->GetYaxis()->GetBinWidth(j)/2 << separator << hist->GetBinError(i,j) << endl;
    if (j>1)
    file_out << endl; // produce a blank line after each set of Ybins for a certain Xbin. Gnuplot likes this.
    }
    if (i>1)
    success=kTRUE;
    }
    file_out.close();
    if (success == kTRUE)
    cout << "*** TRemiHistExport: Histogram " << hist->GetName() << " written to " << filename << endl;
    return success;
    }
    [/code]

I’ve tried to rebin the histogram first and then change the axis’ range but I get the same behaviour.
I’ve also tried

Rebin(4,0)

but again no change.
Any ideas on why is this happening and/or how can it be solved?

2

Either:
TH1F *histo = new TH1F(“hist”, “histo”, channels, 0, channels);
or:
TH1F *histo = new TH1F(“hist”, “histo”, channels, 0.5, (channels + 0.5));
or:
TH1F *histo = new TH1F(“hist”, “histo”, channels, 1, (channels + 1));
and then (always):
histo->Fill(int(x + 0.5), y);

BTW. Note that calling TH1:SetBins(Int_t nx, Double_t xmin, Double_t xmax) will destroy the previous bin contents.

3

Thank you very much for your help.
It doesn’t seem to work however. I tried all the ways

[code]TH1F *histo = new TH1F(“hist”, “histo”, channels, 0, channels);

TH1F *histo = new TH1F(“hist”, “histo”, channels, 0.5, (channels + 0.5));

TH1F *histo = new TH1F(“hist”, “histo”, channels, 1, (channels + 1));[/code]

and then in the filling process

histo->Fill(int(x + 0.5), y);

but it doesn’t seem to change anything. The histogram after rebinning but without scaling the axis is the following

while after x scaling is this one

What do you mean by

Maybe I miss something important. To rebin a histogram “by a factor of 4” means to add all the bin contents every 4 bins and put them in one bin, isn’t that right?

It strikes me strange the fact that, after the rebinning I would expect the final histogram to have bin/4 bins. So in my case the 1st histogram has 4096, which means I would expect the new histogram to have 4096/4=1024 bins, but that doesn’t happen “automatically”. That’s why I am using

TH1:SetBins(Int_t nx, Double_t xmin, Double_t xmax)
4

Attach your data file so that I can try your macro myself.

5

Thank you very much for your time!
Please find attached the input file.
FS3000120Run.dat (132 KB)

6

2048.0000 870.33331
2048.0000 906.00000

7

That’s really weird…
My file comes from the attached fortran code.
I don’t see why this is happening…

REAL E(5000),R(5000),B,C,D,F,H INTEGER Y(5000) OPEN(UNIT=15,FILE='FS3000120Run.dat',STATUS='OLD') DO 40, I=1,4095 READ(15,*) E(I),Y(I) R(I)=Y(I)*1. 40 ENDDO CLOSE(15) OPEN(UNIT=25,FILE='FSS3000120.dat',STATUS='UNKNOWN') DO 60 I=1,288 WRITE(25,*) E(I),R(I) 60 ENDDO DO 50 K=288,700 J=K B=1.30 C=0.70 D=R(K+1)/R(K) IF(D.GE.B.OR.D.LE.C)THEN C F=(R(K-6)+R(K-5)+R(K-4))/3. C H=(R(K+6)+R(K+5)+R(K+4))/3. C R(J+1)=(F+H)/2 R(J+1)=(R(K-2)+R(K-1)+R(K))/3. ENDIF WRITE(25,*) E(J+1),R(J+1) 50 ENDDO DO 70 I=702,4095 WRITE(25,*) E(I),R(I) 70 ENDDO CLOSE(25) STOP END

Also, if I try to rebin it using awk

the result doesn’t have this spike… Now I am really confused…

8

This piece of FORTRAN code does NOT create the “FS3000120Run.dat” file - it reads it and “copies” the data into the “FSS3000120.dat” file (possibly modifying the contents of lines 289 to 701, the remaining ones are simply “copied”).

BTW. The “FS3000120Run.dat” file has two columns with “floating point” numbers, but this piece of FORTRAN code expects the numbers in the second column to be “integer”.

P.S. This piece of AWK code doesn’t give a damn about “$1”, so it’s NOT sensitive to any bugs in “channel numbers”.

9

I tried the same rebinning process(including your corrections) with the attached file and the spike is still there

1280×800

FS3000120Run.dat (132 KB)

10

c1.gif696×502 8.57 KB