Filter and counts events from a TTree

faca87 · December 22, 2022, 6:43pm

Hello, I need a macro to filter events and write the counts in a txt file.
In particular the macro should:

Open the Sim.root file
Read how much events have Edep>0(Edep is a branch and it is a vector). Write the total number of Edep>0 in the txt file
3.Read how much events have Edep[0]>0 (As I said, Edep is a branch and it is a vector. Edep[0] is the first detector). Write the total number of Edep[0]>0 in the txt file
4.Read how much events have Edep[1]>0 (As I said, Edep is a branch and it is a vector. Edep[1] is the second detector). Write the total number of Edep[1]>0 in the txt file
Read how much events have Ep>0 (Ep is other branch) and write the total number of Ep>0 in txt file
Read the total gamma (i.e. the entries of Eg[0] (Eg is a branch and it is a vector).
Read how much Events have Eg[0]>value1 && Eg[0]<value 2, Eg[0]>value3 && Eg[0]<value 3 etc and write these counts in the txt file (Eg is a branch. It is a vectore, but I need Eg[0]).
Copy the Tree in a new file just for Edep>value entries.

I wrote the macro looking this example Root file filter and this tutorial ROOT: tutorials/tree/copytree3.C File Reference but it crashes.

As you can see:

Here I want to count events Edep>0, Edep[0]>0, Edep[1]>0,

  if (Edepentry > 0){
	         DetEvent++;
  		}
  		if (Edep0entry > 0){
	         DetEventPD++;
  		}
  			if (Edep1entry > 0){
	         DetEvent1D++;
  		}

Here I want to fill the TTree when Edep>1000

if (Edepentry > 0){
	         DetEvent++;
  		}

Here I want to counts the total gamma

const auto Eg0branch = tree->GetBranch("Eg[0]");
	 double Eg0total = Eg0branch->GetEntries();

Here I count each gamma of a given energy range

if (Eg0entry > Eg0aMin && Eg0entry <Eg0aMax ){
	        Eg0acount++;
  		}
  		 else if (Eg0entry > Eg0bMin && Eg0entry <Eg0bMax ){
	        Eg0bcount++;
  		}
  		else if (Eg0entry > Eg0cMin && Eg0entry <Eg0cMax ){
	        Eg0ccount++;
  		}
  		else	if (Eg0entry > Eg0dMin && Eg0entry <Eg0dMax ){
	        Eg0dcount++;
  		}
  		else if (Eg0entry > Eg0eMin && Eg0entry <Eg0eMax ){
	        Eg0ecount++;
  		}
  		else if (Eg0entry > Eg0fMin && Eg0entry <Eg0fMax ){
	        Eg0fcount++;
  		}
  		else if (Eg0entry > Eg0gMin && Eg0entry <Eg0gMax ){
	        Eg0gcount++;
  		}

Here I fill the Tree in a new root file when Edep>1000

 if (Edepentry > 1000){
	         copytree->Fill();
  		}

Here macro e root file
treefilter.cpp (3.8 KB)

Wile_E_Coyote · December 22, 2022, 7:09pm

Run:
root -l -q -b Sim.root -e 'Tree1->Print();'
and you will see that all branches are “vector<double>”.

Run:
root -l -q -b Sim.root -e 'Tree1->MakeClass();'
to generate a simple “analysis skeleton”.

Another possibility would be to use the ROOT::RDataFrame (e.g., with the help of the Data Frame tutorials).

faca87 · December 23, 2022, 6:06am

Hi @Wile_E_Coyote

Launching the command on my laboratory machine, I get two files (Tree1.h and Tree1.C). What should I do now to get the information that I need (counts of the branches as written in first post) and to get the new TTree for Edep>1000 ?

Tree1.h (4.3 KB)
Tree1.C (1.4 KB)

BTW: @bellenote I got the two files by laboratory machine, but on my windows I get

c:\my>root -l -q -b Sim.root -e 'Tree1->MakeClass();'
ROOT_cli_1:1:1: warning: multi-character character constant [-Wmultichar]
'Tree1-;'
^
ROOT_cli_1:1:1: warning: character constant too long for its type

c:\my>

Wile_E_Coyote · December 23, 2022, 8:27am

Inside of the “Loop” method, you can use any available std::vector method, e.g.:

unsigned long Edep_size = Edep->size();
double Edep_0 = (Edep_size > 0 ? Edep->at(0) : 0.);
double Edep_1 = (Edep_size > 1 ? Edep->at(1) : 0.);

See also the TTree::MakeClass method description and “ROOT User’s Guide” → “Trees” → “Using TTree::MakeClass”

BTW. If you have problems with single quotation marks, try the double ones:
root -l -q -b Sim.root -e "Tree1->MakeClass();"

faca87 · December 23, 2022, 1:46pm

Thank you @Wile_E_Coyote
but maybe I didn’t understand your suggestion.
That’s because I wrote

To count how many events have been detected by PD detector (i.e. I must count the number of entries Edep[0]>0 )


   double DetEventPD = (Edep_size > 0 ? Edep->at(0) : 0.);

and similarly for the detected events by DD detector i.e. I must count the number of entries Edep[1]>0 )

double DetEventDD = (Edep_size > 0 ? Edep->at(1) : 0.);

To count the number of total emitted primary gammas (i.e. I must counts the total number of entries of Egamma[0])

double PrimGamma = (Egammacascade_size ? Egammacascade->at(0) : 0.);

To count the fusions events (i.e. I must count the number of entries Ep>0)

double FusionEvents = (Ep_size > 0 ? Ep->at(0) : 0.);

To count the number of 1824keV gamma (i.e. I must counts the number of entries with energy 1820<Egammacascade[0]<1828)

double Gamma1824 = (Egammacascade_size > 1820 && Egammacascade_size < 1828 ? Egammacascade->at(0) : 0.);

and similarly for the other energies of gammas (4028keV, 5014keV, etc.)

but as you can see in the txt file, the results are not correct, because, if I look the histograms, I know that, for the small file that I uploaded yesterday:

a. Edep[0]>0 , 1894 entries
b. Edep[1]>0 , 679 entries
c. 1820<Egammacascade[0]<1828 , 226 entries
d. Ep>0 , 23665 entries

Tree1.C (3.5 KB)

Tree1.h (4.3 KB)

Results.txt (311 Bytes)

Wile_E_Coyote · December 23, 2022, 2:11pm

The “(Edep_size > 0 ? Edep->at(1) : 0.)” is wrong.

If you want to access the vector’s element “n”, you must make sure that the vector’s size is at least “n + 1”, e.g.:
(Edep_size > n ? Edep->at(n) : 0.)

See the std::vector::size and the std::vector::at methods descriptions (and try the examples shown therein).

faca87 · December 23, 2022, 2:20pm

Thank you @Wile_E_Coyote
I fixed the error

double DetEventDD = (Edep_size > 1 ? Edep->at(1) : 0.);

but I still get wrong results. Moreover, then I don’t understand how to counts events satisfying some condition…as for example the requests

 Edep[0]>0 
 1820<Egammacascade[0]<1828

Wile_E_Coyote · December 23, 2022, 2:29pm

if ((Edep->size() > 0) && (Edep->at(0) > 0.)) { std::cout << "Edep is fine.\n"; }
if ((Egammacascade->size() > 0) && (Egammacascade->at(0) > 1820.) && (Egammacascade->at(0) < 1828.)) { std::cout << "Egammacascade is fine.\n"; }

faca87 · December 23, 2022, 2:42pm

Sorry @Wile_E_Coyote , probably I didn’t explain what I need (or more probably I’m not understanding your suggestion).

Let’s do an example for Egammacascade[0].

This is the plot of Egammacascade[0] and, as you can see in the statbox, the total number of entries (i.e. the total number of photons) is 23664

Now how want to count how many gamma with energy 1824keV, how many gamma with energy 4028keV etc there are.

For example, for the 1824keV photons, using the TCuts

TString henecut5 = TString::Format("Egammacascade[0]>> htemp(7501, -1., 7500.)");
			TCut cutedep5min = TString::Format("Egammacascade[0] > 1820").Data();;
			TCut cutedep5max = TString::Format("Egammacascade[0] < 1828").Data();;
			te->Draw(henecut5, cutedep5min && cutedep5max );

in the statbox I read 226 entries (i.e. 226 photons have energy 1824keV)

Then, the code that you are helping me to write, should count the number of 1824keV photons (i.e. 226).

Similarly it should count the number of entries for which Edep[0]>0, etc

Wile_E_Coyote · December 23, 2022, 2:53pm

So you can simply define some “counters” (e.g., integer variables) and increment them when the corresponding conditions are met.

faca87 · December 23, 2022, 3:15pm

I’m trying in this way (example to get 1824 keV photons)

for (Long64_t jentry=0; jentry<nentries;jentry++) {
	  	fChain->GetEntry(jentry);
	  	if (Egammacascade->at(0) > 1820 && Egammacascade->at(0) < 1828){
	  		Gamma1824++;
		  }
   }

but when I run t.Loop(); the macro crashes.

@Wile_E_Coyote could you make an example for one of the needed quantity please?
Tree1.C (2.5 KB)

Wile_E_Coyote · December 23, 2022, 4:41pm

What does the “crash message” say?

faca87 · December 23, 2022, 5:00pm

Hi @Wile_E_Coyote …

There isn’t an error message…just it closes ROOT

root [0] .L Tree1.c
root [1] Tree1 t
(Tree1 &) @0x38e0690
root [2] t.Loop();

c:\my>

EDIT I tried by laboratory machine (linux). ON linux I get this error

root [0] .L Tree1.C
root [1] Tree1 t
(Tree1 &) @0x7fdb8738b010
root [2] t.Loop();
Error in <TRint::HandleTermInput()>: std::out_of_range caught: vector::_M_range_check: __n (which is 0) >= this->size() (which is 0)
root [3]

Wile_E_Coyote · December 23, 2022, 5:04pm

So, the message is clear, isn’t it.
You try to access the vector’s element “0”, but the vector’s size is “0”.

faca87 · December 23, 2022, 5:09pm

Sorry @Wile_E_Coyote but I don’t understand one thing. I know that Egammacascade[0] exists, because when I plot, I can plot Egammacascade[0] and, as I shown in my previous message, I could make the Tcuts on it!

Therefore, how can I get the number entries of Egammacascade[0] and the number satisfying some cuts?

Wile_E_Coyote · December 23, 2022, 5:10pm

Use the “full” conditions (checking the vector’s size, too), e.g., like here:

faca87 · December 23, 2022, 5:13pm

I’m sorry for all these questions @Wile_E_Coyote , but in this code

I don’t see a counter! It doesn’t give me the number of entries (i.e… I don’t get 226 photons having energy 1824 keV…)

How can I get the entries?

EDIT I mean that my counter should do just what the statbox does (i.e. counting the entries)

faca87 · December 24, 2022, 10:57am

Hello, I solved myself!

system · January 7, 2023, 10:57am

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.