Make an NTuple for an ADC output file

atha.stam · January 28, 2014, 8:32am

I am trying to sort a bit, the ADC values of a read out system. The output files have the following format

Kmax Event File - Text Format 1 4 1000 65 4121 9426 12312 56 4118 8882 12307 1273 4188 8217 12309 1291 4204 8233 12308 1329 4170 8225 12303 1341 4135 8207 12306 63 4108 8904 12300 60 4106 8897 12307 731 4108 8192 12306 ... ÿÿÿÿÿÿÿÿ

Let me explain more about the file. The 1st line is rather unimportant. The 2nd line has 3 numbers. The first one says the event type(nothing important), the 2nd one is how many detectors the system has(it’s actually the number of the columns and it’s not always 4) while the 3rd one gives the number of events, after which the second line appears.
In this case after 1000 events the line 1 4 1000 will repeat itself. The last line is to be ignored.

The ADC has a number of channels. In the above case it’s 4096(it’s not always like that). So the real number of the 2nd column is “Value of Second Column”-“Number of Channels”, of the 3rd is “Value of Third Column”-2*“Number of Channels” and so on…

Like I said each column represents a detector. My target is to have an NTuple that will sort every detector in one histogram. I need the NTuple, because it will give me the advantage to plot some events, depending on some restrictions/selection over the other detectors.

How can this be achieved? Thank you very much in advance!

Danilo · January 28, 2014, 9:50am

Hi,

this problem is to be solved in my opinion in steps:

Setup the parsing of the file
Filling an ntuple and persist it in a rootfile
Can be achieved coding the rules you describe in C++ using an ifstream (cplusplus.com/reference/fstream/ifstream/). Alternatively, and this is going to be much simpler, using Python.
Once you have in hands a parser for your ADC output the best class you can use for this simple date is the TNtuple. The documentation online should be enough to get you started, as your use case seems to fit exactly the API offered by the class: root.cern.ch/root/html/TNtuple.html

Cheers,
Danilo

atha.stam · January 28, 2014, 10:17am

It could be a good idea to modify the output outside root. I think I can easily do it outside root using awk.
The thing is that it will be much convenient to using on different systems that have only root.

Danilo · January 28, 2014, 10:27am

Hi,

I meant step 1) and 2) to be part of the same application, either written entirely in C++ or in Python using the ROOT bindings, PyRoot.

Cheers,
Danilo

atha.stam · January 28, 2014, 11:38am

That would be a problem, because I don’t know Python…

Danilo · January 28, 2014, 1:37pm

As I said, Python is just an option, the same can be done in C++ easily using the std::ifstream class.

Cheers,
Danilo

atha.stam · January 28, 2014, 1:48pm

I am not sure if I am able to do these stuff in c++.
Once I tried to extract and modify data from a relevant file and it was rather difficult.
Back then I discovered awk, which made my life much easier.
I suppose there isn’t any relevant option on root,is there?

atha.stam · January 31, 2014, 3:32pm

I found a way to manipulate my data. I am using this very script

[code]#!/bin/bash

if test $1 ; then
if [ -f $1.evnt ] ; then
rm -f $1.dat
sed -n ‘2p’ $1.evnt | (read v1 v2 v3
for filename in $1*.evnt ; do
echo -e “Processing file $filename"
sed ‘$d’ < $filename > $1_tmp
sed -i ‘/Kmax/d’ $1_tmp
sed -i ‘/^’”$v1"’ ‘"$v2"’ /d’ $1_tmp
cat $1_tmp >> $1.dat
done
v3=wc -l $1.dat | awk '{print $1}'
echo -e “$v1 $v2 $v3” > .$1.dat
rm -f $1_tmp)
else
echo -e "\a!!!"
echo -e " Event file $1.evnt doesn’t exist !!!"
echo -e "!!!"
fi
else
echo -e "\a!!!"
echo -e "!!! Give name for event files !!!"
echo -e "!!!"
fi
awk ‘{print $1, $2-4096, $3-(24096), $4-(34096)}’ $1.dat >$1_Processed.dat
rm -f $1.dat
exit 0[/code]

So I used root’s tutorial to make the ntuple and some histograms. The macro is the following

[code]#include "Riostream.h"
void ntuple() {

gSystem->Exec("./evnt2dat c+Au");
TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
dir.ReplaceAll(“ntuple.C”,"");
dir.ReplaceAll("/./","/");
ifstream in;
in.open(Form("%sc+Au_Processed.dat",dir.Data()));

Float_t de1,e1,de2,e2;
Int_t events = 0;
TFile *f = new TFile(“ntuple.root”,“RECREATE”);
//Create the histograms
TH1F *histo_de1 = new TH1F(“histo_de1”,“x distribution”,4096,1,4096);
TH1F *histo_e1 = new TH1F(“histo_e1”,“x distribution”,4096,1,4096);
TH1F *histo_de2 = new TH1F(“histo_de2”,“x distribution”,4096,1,4096);
TH1F *histo_e2 = new TH1F(“histo_e2”,“x distribution”,4096,1,4096);
TH2F *dee = new TH2F(“dee”, “de VS e”,4096,1,4096,4096,1,4096);
//Create the ntuple
TNtuple *ntuple = new TNtuple(“ntuple”,“data from ascii file”,“de1:e1:de2:e2”);

//Fill the histos and the ntuple
while (1) {
in >> de1 >> e1 >> de2 >> e2;
if (!in.good()) break;
//if (nlines < 5) printf(“x=%8f, y=%8f, z=%8f\n”,x,y,z);
histo_de1->Fill(de1);
histo_e1->Fill(e1);
histo_de2->Fill(de2);
histo_e2->Fill(e2);
dee->Fill(de1, e1);
ntuple->Fill(de1,e1,de2,e2);
events++;
}
printf(" found %d events\n",events/4);

in.close();

f->Write();
//h1->Draw();
dee->Draw(“COLZ”);
}[/code]

You can see that using this way, I will have to give root every time I execute the macro, the name of the file to process. Is there a way to do this automatically?
I also found that the COLZ option makes the built last a lot longer. Is this typical? Is there a way to make it faster?

atha.stam · February 1, 2014, 3:10pm

I was thinking of something like

Can something like this this be achieved?

jfcaron · February 2, 2014, 7:08am

If you “.L” your macro to get a function that you can call at the ROOT prompt, then you can give it arguments. If you make your function take a TString or char * argument, then you can give it the name of a file in double-quotes, like a normal function.

From the command line (outside of root), if your macro can take arguments then you can do:

Note the quotes around the macro name and its brackets. If you wish to pass a string to the macro this way, it needs to be in double-quotes too, but the macro name and arguments need to be in quotes, so you might need to do some escaping with /.

Take a look here: http://bazaar.launchpad.net/~jfcaron/+junk/TRIUMFBeamTest/view/head:/standard_analysis/standard_analysis.C at lines 2-3 and 62-67 for some example code that I used.

Jean-François

atha.stam · February 2, 2014, 9:23am

Thank you very much for your help.
The thing is that I am not sure, if I made clear what I want.
My macro includes the following commands

gSystem->Exec("./evnt2dat c+Au"); in.open(Form("%sc+Au_Processed.dat",dir.Data()));

The c+Au is the file the script evnt2dat uses to produce the c+Au_Processed. The last file, is the ascii file that builds the ntuple.

So I was thinking something like

$ root ntuple.C

where ntuple.C

... gSystem->Exec("./evnt2dat <filename>"); ... in.open(Form("%s<filename>_Processed.dat",dir.Data())); ...

jfcaron · February 2, 2014, 7:18pm

You’ll have to do something like:

and your ntuple.C file looks like this:

void ntuple(char * fname_c)
{
TString fname(fname_c);
...
gSystem->Exec(TString::Format("./evnt2dat %s", fname.Data());
...
in.open(TString::Format("%s%s_Processed.dat",dir.Data(),fname.Data()));
...
}

Jean-François

atha.stam · February 2, 2014, 8:48pm

Thank you very much for your help!!!
It really works!

I have some questions though, on how how it works…

TString::Format : I assume using it, cint expects to read characters.
%s : Tells to read the variable that follows next
%s%s : It is connected to dir.Data() and fname.Data()(first %s for dir and second for fname. What does .Data() do?

The tricky part is the following: As you have noticed I am creating 4(+2) histograms. Each 1D for every column and one 2D per 2 columns. The number 4, comes from the columns of the input file. I am able to create a variable in my script, that will show the number of columns, therefore the number of histos. The thing is, will root be able to create those histos? Can TH1F and TH2F be created inside a for loop?

jfcaron · February 3, 2014, 2:03am

TString::Format is described here: http://root.cern.ch/root/html/TString.html#TString:Format, it uses the same codes as the old printf function: http://www.cplusplus.com/reference/cstdio/printf/.

The %s code means “string”, but it means a C-style string like char *, to the .Data() is necessary to get a “C-style” version of the TString. You could use the fname_c directly, but I always convert to TStrings to make manipulations easier. The parameters are in sequential order, so %s%s tells Format to expect two string arguments. Play around with it if you are unfamiliar.

You can create histograms (or any object) in a loop. You can do something like:

std::vector<TH1F> myhists;
for(Int_t i=0;i<4;i++)
{
  TString name = TString::Format("h_%d",i);
  TH1F temp(name,name,10,0,10);
  myhists.push_back(temp);
}

You need to give the histograms different names because of the way ROOT manages objects.

atha.stam · February 3, 2014, 7:34am

Thank you once again for your info!
There is still a problem though(if you feel that I should create a new topic on how to create histos inside a for loop please let me know).

I used your code-after commenting my histos constructors- but I get the following error

Then I tried to play a bit, creating the histos just how I want them, but something feels wrong. My parts of the code are

... Float_t de1,e1,de2,e2; Int_t channels=4096; ... //Create the e histos std::vector<TH1F> my_e_hists; for(Int_t i=1;i<3;i++) { TString name = TString::Format("histo_e%d",i); TH1F temp(name,name,channels,1,channels); my_e_hists.push_back(temp); } //Create the de histos std::vector<TH1F> my_de_hists; for(Int_t j=1;j<3;j++) { TString name = TString::Format("histo_de%d",i); TH1F temp(name,name,channels,1,channels); my_de_hists.push_back(temp); } ... //Fill the histos and the ntuple while (1) { in >> de1 >> e1 >> de2 >> e2; if (!in.good()) break; for(Int_t k=1;k<3;k++) { histo_e$k->Fill(e_k); } for(Int_t l=1;l<3;l++) { histo_e$l->Fill(de_l); } //No idea how to do a for loop to fill the ntuple... ntuple->Fill(de1,e1,de2,e2); ... }

At this point I feel lost…(and I should open a new thread…)