Efficient loading of many ".dat" files and plotting them as Histograms

myless · November 4, 2019, 12:00am

ROOT Version: 6.18/04
Platform: Arch Linux-64 bit
Compiler: linuxx8664gcc

Hello everyone, I am measuring the Compton effect using a Cs137 Source and NaI detectors. The output of the data consists of multiple “.dat” files with four columns of data. The column I am interested in is the third column where the “channel” is listed. I essentially want to plot a histogram of the frequency of the channels and use this histogram to calibrate my experiment, i.e find out that channel number 12000 refers to 662keV, and channel 2000 refers to 32keV, etc.

I have 20 of these “.dat” files and have previously used this quick/dirty method in RInt:

tr = TTree("tr", "title");
tr.ReadFile("filename","a:b:c:d");
tr.Draw("c>>hist(bins,0,maxchannel)");

What I get is a Histogram, but not only is this tedious to repeat for 20 or so files, but I can’t figure out how to label specific peaks in my data, for calibration purposes. I can however right click on the histogram and select the ShowPeaks option in the GUI Plotter.

My questions are:

Is there a standard way in ROOT to load all of the files in a certain directory into ROOT and plot them as histograms? (In Python I would just do a forloop with plt.hist(), but I am not sure how I would achieve the same effect in ROOT).
Is there a way to label and display the specific bin/channel peak observed in the histogram?

Here is an example of what I am plotting individually now.

Thank you,
Myles

bellenot · November 4, 2019, 8:30am

Maybe something like this:

   TString filename;
   void* dirp = gSystem->OpenDirectory(dirname);
   if (dirp) {
      Int_t n = 0;
      while ((filename = gSystem->GetDirEntry(dirp)) != "" ) {
         if (filename.EndsWith(".dat")) {
            tr = TTree("tr", "title");
            tr.ReadFile("filename","a:b:c:d");
            tr.Draw(Form("c>>hist%d(bins,0,maxchannel)", ++n));
         }
      }
   }

Then you can loop over your histograms and call TH1F::ShowPeaks().

Maybe @couet can give an answer to that question…

couet · November 4, 2019, 8:53am

You can add a TText at the peak position.

   auto l = new TText( x_peak_pos, y, Form("%g",x_peak_pos));
   l->Draw();

myless · November 5, 2019, 1:37am

Thank you both for your answers.

I have modified the code for my specific application:

#include <iostream>
 #include <string>
 
 void fileload(){
     //Variables
     TString dirname = "data/";
     TString filename;
     Int_t bins = 200;
     Int_t maxchannel = 30000;
 
 
 
     //beginning of fileloading
     void* dirp = gSystem->OpenDirectory(dirname);
     if (dirp) {
         Int_t n = 0;
         //TTree tr = TTree("tr","title");
         while ((filename = gSystem->GetDirEntry(dirp)) != ""){
             if (filename.EndsWith(".dat")){
                 TTree tr = TTree("tr","title");
                 tr.ReadFile("filename","a:b:c:d");
                 tr.Draw(Form("c>>hist%d(bins,0,maxchannel)", ++n));
                 
             }
         }
     }   
 }

However, when I run this code in RInt, using:

root -l .x fileload.C

I get the following errors :

Processing fileload.C...
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename
Error in <TTree::ReadFile>: Cannot open file: filename

Good news seems to be that the files are noticed in the directory, bad news is that the files are not being loaded.

As for the Peak Positions, I am unsure what the x_peak_pos and associated y are, I tried:

TSpectrum *s=new TSpectrum(2,1);
int n_peaks=s->Search("hist",0.03,"",0.5);
double* peak_position_array=s->GetPositionX();

Where this can return an array with the x-value of a peak, for example, with the hist file attached, I can replace “hist” with “h__1” and type:

peak_position_array[0]

after the previous lines of code to get the peak x-value.

Is there a more efficient method?I am curious how the GUI is able to quickly find the peaks and mark them. Perhaps there is a way for the GUI to spit out the x and y values of the markers?

Thank you both for your help,

Myles

bellenot · November 5, 2019, 8:26am

Sorry, my bad. Please replace tr.ReadFile("filename","a:b:c:d"); by tr.ReadFile(filename.Data(),"a:b:c:d");

couet · November 5, 2019, 9:22am

myless:

As for the Peak Positions, I am unsure what the x_peak_pos and associated y are, I tried:
TSpectrum *s=new TSpectrum(2,1);
int n_peaks=s->Search("hist",0.03,"",0.5);
double* peak_position_array=s->GetPositionX();

You get an array with the X positions of all the peaks. That’s why it is an array. That was your question ?

myless · November 5, 2019, 2:07pm

Bertrand, Olivier,
Sorry for the poor wording of my question.

I have figured out how to get an array with the X-Values of the peaks, but I cannot figure out how to get a y-value associated to each of them.

From reading this:

And from the documentation of TText: TText

It seems that I not only need the x-values, (which I can find now), but that I also need the y-value associated with the peak. I have attached the code and an example .dat file of how I am finding the peaks. peakfinder.C (537 Bytes)
Link to Data: Link to DataSet Folder
The filename variable may need to be changed to get it to work.
An example of what the code can do on my computer is:

With an output of:

Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
Warning in <TSpectrum::SearchHighRes>: Peak buffer full
The number of peaks found are: 4
12510
810
4190
3110

Thank you for your suggestion. Implementing your change now spits out the correct name of each file, but still results in an error. Here is the error:

Processing fileload.C...
Error in <TTree::ReadFile>: Cannot open file: coin-cs-60_001_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-90_001_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: coincidence-offcenter-Na_004_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-60_001_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: norm-cal-Cs-1.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-15_001_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: norm-cal-Na-1.dat
Error in <TTree::ReadFile>: Cannot open file: coincidence-offcenter-Na_004_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-120_001_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: coincidence-center-Na_004_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-90_001_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-120_001_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-15_001_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: norm-cal-Cs-0.dat
Error in <TTree::ReadFile>: Cannot open file: norm-cal-Na-0.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-30_001_ls_0.dat
Error in <TTree::ReadFile>: Cannot open file: coincidence-center-Na_004_ls_1.dat
Error in <TTree::ReadFile>: Cannot open file: coin-cs-30_001_ls_1.dat

Thank you both,

Myles

couet · November 5, 2019, 2:57pm

once you have the x position xp you can retrieve the Y value with:

      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->GetBinContent(bin);

bellenot · November 5, 2019, 3:18pm

Here is how to properly read your files, fill the histograms, and display the canvas (I tried with two of your files):

#include "TSystem.h"
#include "TString.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TTree.h"
#include <vector>
#include <iostream>
using namespace std;

void fileload()
{
   TString dirname = "DataSet/";
   TString filename;
   Int_t bins = 200;
   Int_t maxchannel = 30000;
   std::vector<TH1F*> hist;
   std::vector<TCanvas*> canvas;

   void* dirp = gSystem->OpenDirectory(dirname);
   if (dirp) {
      Int_t n = 0;
      while ((filename = gSystem->GetDirEntry(dirp)) != "") {
         if (filename.EndsWith(".dat")) {
            filename.Prepend(dirname);
            std::cout << "Reading " << filename.Data() << std::endl;
            hist.push_back(new TH1F(Form("hist[%d]",n), Form("hist[%d]",n), 1000, 0, 20000));
            TTree tr("tr","title");
            tr.ReadFile(filename.Data(), "a:b:c:d");
            canvas.push_back(new TCanvas(Form("canvas%d",n), Form("canvas%d",n)));
            canvas[n]->cd();
            tr.Draw(Form("c>>hist[%d]",n));
            canvas[n]->Modified();
            canvas[n]->Update();
            ++n;
         }
      }
   }
}

myless · November 5, 2019, 3:42pm

Thank you both so much for the help. This code works well; I am shocked how fast it is… I will fiddle around with the syntax to implement the TSpectrum finder as well.

One final question if I may, what is the advantage of including the T functions at the beginning of the program? I find that my macros still compile even without them. If this is easily answered somewhere else, then could I please be pointed to that resource?

Best,

Myles

bellenot · November 5, 2019, 3:45pm

Well, in principle it’s required when compiling it (I use ACLiC), like any standard C++ code

bellenot · November 5, 2019, 3:57pm

BTW, here is another version, with separated loop for the canvas creation and drawing:

#include "TSystem.h"
#include "TString.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TTree.h"
#include <vector>
#include <iostream>
using namespace std;

void fileload()
{
   TString dirname = "DataSet/";
   TString filename;
   Int_t n = 0;
   Int_t bins = 200;
   Int_t maxchannel = 30000;
   std::vector<TH1F*> hist;
   std::vector<TCanvas*> canvas;

   void* dirp = gSystem->OpenDirectory(dirname);
   if (dirp) {
      while ((filename = gSystem->GetDirEntry(dirp)) != "") {
         if (filename.EndsWith(".dat")) {
            filename.Prepend(dirname);
            std::cout << "Reading " << filename.Data() << std::endl;
            hist.push_back(new TH1F(Form("hist[%d]",n), Form("hist[%d]",n), 1000, 0, 20000));
            TTree tr("tr","title");
            tr.ReadFile(filename.Data(), "a:b:c:d");
            tr.Draw(Form("c>>hist[%d]",n));
            ++n;
         }
      }
   }
   /// remove this loop if canvases and drawing are not needed
   for (int i=0;i<n;++i) {
      canvas.push_back(new TCanvas(Form("canvas%d",i), Form("canvas%d",i)));
      hist[i]->Draw();
   }
}

myless · November 6, 2019, 3:38am

Ahh I see. I will give this a try to practice better habits. Thank you for the explanation.