Help on Method/s to create Tree from multi-fold gamma-ray coincidence data

ROOT
47

After building almost all the tools required for multi-fold gamma-ray coincidence data analysis, today I tried to sort ~ 8 GB data (28 data files) into a tree (“canSort”) using “can2root.cxx”. This has generated a ROOT file “sixT2MeV.root” (3.2 GB). When I use this ROOT file later to create a symmetric 3D THnSparse (“his3D”) from “gen_Cub.cxx” with “CubDim = 4096”, I am getting the following warning:

Warning: MaxNbins reached for the his3D.

his3D->Print(“all”); gives

THnSparseT<TArrayF> (*0x19a0590): "his3D" "his3D"
  3 dimensions, 1.71278e+07 entries in 16777217 filled bins
    axis 0 "EclabX": 4096 bins (0..4096), variable bin sizes
    axis 1 "EclabY": 4096 bins (0..4096), variable bin sizes
    axis 2 "EclabZ": 4096 bins (0..4096), variable bin sizes

Does the above warning imply loss of the data/events?
It is not clear to me where and how to find description on limiting values of possible number of bins for THnSparse. I know, in our case it is CubDim * CubDim. But why?

48

Yes, this warning means that some data have been skipped.

I set “MaxNbins = CubDim * CubDim;” just to limit the amount of used RAM to something small (i.e. 64MB in this case) but, you can set whatever you want, e.g. “MaxNbins = 268435454;” (at least around 1GB RAM may be needed when a bin content is held by a Float_t).

49

More methods’ headers are now in the doxygen format.

50

Strange Behavior of different types of Projections:

I am trying to understand behavior and performance of TH3D and THnSparse. I have generated TWO different ROOT file: one using TH3D and another using THnSparse; with the same conditions. I used the 3D histogram (his3D) in them to generate 1D projection using different methods as follows:

==============================================================
METHOD1: (using TH3D Class)

TFile *f = new TFile("cub_his3D.root");
his3D->GetXaxis()->SetRangeUser(668,672);
his3D->GetYaxis()->SetRangeUser(449,453);
his3D->Project3D("z")->Draw();

This gives me 16566 entries in a resulting histogram.

METHOD2: (using TH3D Class)

TFile *f = new TFile("cub_his3D.root")
his3D->SetAxisRange(668,672,"X");
his3D->SetAxisRange(449,453,"Y");
TH1D *pjz1 = his3D->Project3D("Z");

This gives me 19268 entries in a resulting histogram.

METHOD3: (using THnSparse Class)

TFile *f = new TFile("cub_sparse.root");
his3D->GetAxis(0)->SetRangeUser(668,672);
his3D->GetAxis(1)->SetRangeUser(449,453);
TH1D *pjz = his3D->Projection(2);
pjz->Draw();

In this case the resulting 1D-histogram has ONLY 10663 entries.

METHOD4: (using THnSparse Class)

TFile *f = new TFile("cub_sparse.root");
TAxis *a0 = his3D->GetAxis(0);
TAxis *a1 = his3D->GetAxis(1);
a0->SetRange(a0->FindFixBin( 668 ), a0->FindFixBin( 672 ));
a1->SetRange(a1->FindFixBin( 449 ), a1->FindFixBin( 453 ));
TH1D *pjz = his3D->Projection(2);
pjz->Draw();

Now the total entries in the 1D-histogram are 12435.

==============================================================

You see, there is a problem! The entries vary from 10663 to 19268 depending on the method of projection and the Class Type. I get more counts in 2nd method. But the problem (as discussed in earlier posts) that I cannot create a TH3D with large dimensions. For the present exercise, I have limited the dimensions to 620 only.

Is the above behavior expected? Which one of the above methods should I trust?

51

For “METHOD1” and “METHOD2” you should use:

TFile *f = TFile::Open("cub_his3D.root");
TH3 *his3D; f->GetObject("his3D", his3D);

For “METHOD3” and “METHOD4” you should use:

TFile *f = TFile::Open("cub_sparse.root");
THnSparse *his3D; f->GetObject("his3D", his3D);

I’m not sure if the “number of entries” is what you want. I guess you’d better try with the TH1::Integral method.

Note also that, in any operations, entire bins are always taken into account so, if you want to compare numbers, you need to make sure that e.g. zaxis->GetBinLowEdge(zaxis->FindFixBin( zmin )) and zaxis->GetBinUpEdge(zaxis->FindFixBin( zmax )) (and so on) are equal for all histograms that you want to compare (btw. following the standard convention for numbering bins, the low-edge is included but the upper-edge is excluded).

Please compare input parameters of the TAxis::SetRange, the TAxis::SetRangeUser and the TH1::SetAxisRange methods.

52

Hi,
I have tried the methods that you have suggested. But still the problem persists. May be, I am not properly following what you said. Anyway, here is what I do:

======================

TFile *f = TFile::Open("cub_his3D.root");
TH3 *his3D; f->GetObject("his3D", his3D);

his3D->GetXaxis()->SetRangeUser(668,672);
his3D->GetYaxis()->SetRangeUser(449,453);
TH1D* pjz1 = (TH1D*)his3D->Project3D("z");
Int_t xlow = pjz1->GetXaxis()->GetBinLowEdge(pjz1->FindFixBin(100));
Int_t xhi = pjz1->GetXaxis()->GetBinUpEdge(pjz1->FindFixBin(700));
cout << pjz1->Integral(xlow,xhi) << endl;
12210

his3D->SetAxisRange(668,672,"X");
his3D->SetAxisRange(449,453,"Y");
TH1D *pjz2 = (TH1D*)his3D->Project3D("z");
Int_t xlow2 = pjz2->GetXaxis()->GetBinLowEdge(pjz2->FindFixBin(100));
Int_t xhi2 = pjz2->GetXaxis()->GetBinUpEdge(pjz2->FindFixBin(700));
cout << pjz2->Integral(xlow2,xhi2) << endl;
14186

======================

53

Just to exclude floating point “slips”, could you, please, try with “668.1, 671.9”, “449.1, 452.9”, “100.1” and “699.9” (btw. are you sure these ranges do not go into underflow or overflow bins?).

54

Wow!! Your last suggestion worked. Thank you.

Here are the results:

===========================

root [0] TFile *f = TFile::Open("cub_his3D_old.root");
root [1] TH3 *his3D; f->GetObject("his3D", his3D);
root [2] his3D->GetXaxis()->SetRangeUser(668.1,671.9);
root [3] his3D->GetYaxis()->SetRangeUser(449.1,452.9);
root [4] TH1D* pjz1 = (TH1D*)his3D->Project3D("z");
root [5] Int_t xlow = pjz1->GetXaxis()->GetBinLowEdge(pjz1->FindFixBin(100.1));
root [6] Int_t xhi = pjz1->GetXaxis()->GetBinUpEdge(pjz1->FindFixBin(699.9));
root [7] cout << pjz1->Integral(xlow,xhi) << endl;
12210

root [8] his3D->SetAxisRange(668.1,671.9,"X");
root [9] his3D->SetAxisRange(449.1,452.9,"Y");
root [10] TH1D *pjz2 = (TH1D*)his3D->Project3D("z");
root [11] Int_t xlow2 = pjz2->GetXaxis()->GetBinLowEdge(pjz2->FindFixBin(100.1));
root [12] Int_t xhi2 = pjz2->GetXaxis()->GetBinUpEdge(pjz2->FindFixBin(669.9));
root [13] cout << pjz2->Integral(xlow2,xhi2) << endl;
12210

===========================

So which one should I believe?

It would be very cumbersome for us to use “668.1” instead of “668” and so on all the time, given that we have to generate hundreds of spectra with different gamma ray windows.

55

Well, you don’t need ±0.1. A small value like ±1.e-6 (or even ±1.e-9) would be sufficient (but, knowing that the ADC channel’s width is 1, I would still prefer ±0.5, of course, so I would use “668.5, 671.5”, “449.5, 452.5”, “100.5” and “699.5”).

BTW. @moneta If you look into the TAxis::SetRangeUser source code, you will find some “fixes for numerical error and for https://savannah.cern.ch/bugs/index.php?99777” inside, which are not present in the TH1::SetAxisRange source code. Maybe one should make both methods “consistent”?

56

What if I add “0.5” to “x, y, z” before filling the 3D histogram, as suggested by you earlier? Will I still need to use 668.5 etc. at a later stage while generating 1D histograms as discussed in earlier posts?

57

Correcting raw ADC values is one thing and bin edges’ “slips” is another one (but they are both related to floating point “rounding errors”, of course). In BOTH cases, if you want to be sure that everything is kosher, you need at least a small shift, like ±1.e-6 (or even ±1.e-9). Then, if you want, you can go to ±0.5 but, you should discuss the +0.5 raw ADC value correction with your colleagues as they maybe do not use it and then you may possibly get some results shifted by +0.5 as compared to theirs (a +1.e-6 or +1.e-9 shift would be completely negligible, of course).

58

Assuming the methods in the post No. 48 of this discussion are correct, we can conclude that we should get 12210 counts between bin no. 100 & 700. While checking if there are any other simpler ways of doing the same, I found the following methods, which employ TAxis::SetRangeUser are also correct. Here we don’t have to bother about floating points corrections, and “integers” can be used to set the ranges.

=============================

TFile *f1 = TFile::Open("cub_his3D_old.root");
TH3 *his3D; f1->GetObject("his3D", his3D);
his3D->GetXaxis()->SetRangeUser(668,672);
his3D->GetYaxis()->SetRangeUser(449,453);
TH1D* pjz1 = (TH1D*)his3D->Project3D("z");
cout << pjz1->Integral(100,700) << endl;
12210

However, there is ONE PROBLEM, the 1D-histogram (“pjz1”) created above is NOT listed in the histogram list. Please see below.

.ls
TFile**         cub_his3D_old.root
 TFile*         cub_his3D_old.root
  OBJ: TH3F     his3D   his3D : 0 at: 0x9f92850
  OBJ: TH1D     his3D_z his3D z projection : 0 at: 0x9d4b6e0
  KEY: TH3F     his3D;1 his3D

=============================

TFile *f = TFile::Open("cub_sparse_old.root");
THnSparse *sparse3D; f->GetObject("his3D", sparse3D);
sparse3D->GetAxis(0)->SetRangeUser(668,672);
sparse3D->GetAxis(1)->SetRangeUser(449,453);
TH1D *pjz = sparse3D->Projection(2);
cout << pjz->Integral(100,700) << endl;
12210

Again, the same PROBLEM is observed, the 1D-histogram (“pjz”) created above is NOT listed in the histogram list. Please see below.

.ls
TFile**         cub_sparse_old.root
 TFile*         cub_sparse_old.root
  OBJ: TH1D     his3D_proj_2    his3D projection EclabZ : 0 at: 0xa82e0c0
  KEY: THnSparseT<TArrayF>      his3D;1 his3D

=============================

What are the methods to “see” pjz and pjz1 which are 1D projections of 3D histogram and THnSparse?

59

No idea what you’re talking about.
I can clearly see “his3D_z” and “his3D_proj_2”. Try:

std::cout << pjz->GetName(); << std::endl;

BTW. It seems to me that you use wrong parameters when calling the TH1::Integral method (I missed this problem in your recent posts). You should probably use:

pjz->Integral(pjz->FindFixBin(100.5), pjz->FindFixBin(699.5))
60

I am sorry, may be I didn’t make my point clear. Let me try again.

(1) In the first instance, I do the following:

root [0] TFile *f1 = TFile::Open("cub_his3D_old.root");
root [1] TH3 *his3D; f1->GetObject("his3D", his3D);
root [2] his3D->GetXaxis()->SetRangeUser(668,672);
root [3] his3D->GetYaxis()->SetRangeUser(449,453);
root [4] TH1D* pjz1 = (TH1D*)his3D->Project3D("z");
root [5] .ls
TFile**		cub_his3D.root	
 TFile*		cub_his3D.root	
  OBJ: TH3F	his3D	his3D : 0 at: 0x1f8c450
  OBJ: TH1D	his3D_z	his3D z projection : 0 at: 0x1f8d8b0
  KEY: TH3F	his3D;1	his3D
root [6] pjz1->Draw();

“pjz1” is a spectrum of transitions which are in coincidence with 670 & 451 keV gamma transitions.

For further analysis. I have to check several such coincidences. For example, let us generate a coincidence spectrum between 670 & 202 as follows:

root [7] his3D->GetYaxis()->SetRangeUser(200,204);
root [8] TH1D* pjz2 = (TH1D*)his3D->Project3D("z");
root [9] .ls
TFile**		cub_his3D.root	
 TFile*		cub_his3D.root	
  OBJ: TH3F	his3D	his3D : 0 at: 0x1f8c450
  OBJ: TH1D	his3D_z	his3D z projection : 0 at: 0x1f8d8b0
  KEY: TH3F	his3D;1	his3D
root [10] pjz1->Draw();
root [11] pjz2->Draw();

I can see “his3D_z” which is same as “pjz2”. But when I display “pjz1”, it is EXACTLY same as “pjz2”.
I have LOST my first spectrum i.e. “pjz1”.

Do you now see my problem? If I generate several such projections, I will always be left with the spectrum corresponding to the last projection. All my earlier spectra are overwritten by the last one.

61

I checked. It seems, both the methods yield the same results.

root [0] TFile *f1 = TFile::Open("cub_his3D.root");
root [1] TH3 *his3D; f1->GetObject("his3D", his3D);
root [2] his3D->GetXaxis()->SetRangeUser(668,672);
root [3] his3D->GetYaxis()->SetRangeUser(449,453);
root [4] TH1D* pjz1 = (TH1D*)his3D->Project3D("z");
root [5] cout << pjz1->Integral(100,700) << endl;
10854
root [6] Int_t xlow = pjz1->GetXaxis()->GetBinLowEdge(pjz1->FindFixBin(100));
root [7] Int_t xhi = pjz1->GetXaxis()->GetBinUpEdge(pjz1->FindFixBin(700));
root [8] cout << pjz1->Integral(xlow,xhi) << endl;
10854
root [9]

Please note that, I am using different 3D-histogram (“cub_his3D.root”) since I don’t have access to the old one (“cub_his3D_old.root”) at present. Hence, slightly different counts.

62

You need to manually rename your histograms right after they are created. Something like:

pjz1->SetName("his3D_670_451");
// ...
pjz2->SetName("his3D_670_202");

BTW. You still use wrong TH1::Integral parameters (see my previous post again).

1 Like
63

Thank you. And yes, I understand my mistake while setting parameters of the TH1::Integral function.

64

After successful testing of the macro on a smaller data (Tree of size 2.3 GB) to generate a THnSparse, I tried to run it on the full data set (~ 18 GB). But it didn’t run successfully. I am attaching the macro as well as the output at the end of execution below.

root [0] .x gen_Sparse3D.cxx++
Info in <TUnixSystem::ACLiC>: creating shared library /home/ajay/11B208Pb_Jul2016/ROOT_Sort/Sparse3D/./gen_Sparse3D_cxx.so
THnSparseT<TArrayF> (*0x2ea0030): "his3D" "his3D"
  3 dimensions, 1.64741e+08 entries in 128708567 filled bins
    axis 0 "EclabX": 2048 bins (0..2048), variable bin sizes
    axis 1 "EclabY": 2048 bins (0..2048), variable bin sizes
    axis 2 "EclabZ": 2048 bins (0..2048), variable bin sizes
0 : { 940.5 , 697.5 , 71.5 } : 1
1 : { 940.5 , 71.5 , 697.5 } : 1
2 : { 697.5 , 940.5 , 71.5 } : 1
3 : { 697.5 , 71.5 , 940.5 } : 1
4 : { 71.5 , 940.5 , 697.5 } : 1
5 : { 71.5 , 697.5 , 940.5 } : 1
6 : { 2048.5 , 1286.5 , 165.5 } : 1
7 : { 2048.5 , 165.5 , 1286.5 } : 1
8 : { 1286.5 , 2048.5 , 165.5 } : 1
9 : { 1286.5 , 165.5 , 2048.5 } : 1
Error in <TBufferFile::WriteByteCount>: bytecount too large (more than 1073741822)
Error in <TBufferFile::WriteByteCount>: bytecount too large (more than 1073741822)
Error in <TBufferFile::WriteByteCount>: bytecount too large (more than 1073741822)

You will notice that I have already reduced the THnSparse dimensions to 2048 from 4096.

I can generate separate THnSparse for smaller data sets, but then there is no provision to add THnSparse, I guess.

Please suggest - how to proceed in this situation.

BTW, I am using 64-bit Fedora20, RAM: 32GB, 24 Cores , ROOT 6.15/01 (linuxx8664gcc)

gen_Sparse3D.cxx (5.7 KB)

65

The “his3D” is a “THnSparseF” and you do not use “weights” when filling it so, it should be 4 bytes per bin.
There are 128708567 filled bins so, it should be around 491MB (which is significantly less than 1GB).
Why does it need more than 1 GB then?
Maybe the chunks’ management uses a lot of memory, too?
I’m afraid we need @pcanal for further debugging.

@pcanal … Can one somehow write this histogram in “split” mode (so that the 1GB limit is valid for “independent parts” / “data members” of it only)?

In the meantime, you may try the following “brutal fix” …

  // TFile *f = TFile::Open("cub_sparse_fifT7.root", "RECREATE", "", 101);
  f->cd();
#if 0 /* 0 or 1 */
  his3D->Write(); // warning: we need to "manually" save the THnSparse
#else /* 0 or 1 */
  // save the "his3D" in a tree ...
  TTree *t = new TTree("tree3D", "tree with his3D");
  t->Branch("his3D.", his3D, 32000, 111); // "max" splitting?
  t->Fill();
  t->Write();
  delete t;
#endif /* 0 or 1 */
  f->Write();
  delete f;

and if it goes fine, you can retrieve the histogram using:

{
  TFile *f = TFile::Open("cub_sparse_fifT7.root");
  TTree *t; f->GetObject("tree3D", t);
  THnSparse *his3D = 0;
  t->SetBranchAddress("his3D.", &his3D);
  t->GetEntry(0);
  delete f; // automatically deletes "t", too
  his3D->Print("all"); // "his3D" stays in RAM
}
66

Hi,

I tried your suggestion. But it failed! However, this time Error message was displayed only once instead of 3 times (if that means anything for debugging) at the end. Please see below:

root [2] .x gen_Sparse3D.cxx++
Info in <TUnixSystem::ACLiC>: creating shared library /home/ajay/11B208Pb_Jul2016/ROOT_Sort/Sparse3D/./gen_Sparse3D_cxx.so
THnSparseT<TArrayF> (*0x3fb8740): "his3D" "his3D"
  3 dimensions, 1.64741e+08 entries in 136400687 filled bins
    axis 0 "EclabX": 4096 bins (0..4096), variable bin sizes
    axis 1 "EclabY": 4096 bins (0..4096), variable bin sizes
    axis 2 "EclabZ": 4096 bins (0..4096), variable bin sizes
0 : { 940.5 , 697.5 , 71.5 } : 1
1 : { 940.5 , 71.5 , 697.5 } : 1
2 : { 697.5 , 940.5 , 71.5 } : 1
3 : { 697.5 , 71.5 , 940.5 } : 1
4 : { 71.5 , 940.5 , 697.5 } : 1
5 : { 71.5 , 697.5 , 940.5 } : 1
6 : { 2105.5 , 1286.5 , 165.5 } : 1
7 : { 2105.5 , 165.5 , 1286.5 } : 1
8 : { 1286.5 , 2105.5 , 165.5 } : 1
9 : { 1286.5 , 165.5 , 2105.5 } : 1
Error in <TBufferFile::WriteByteCount>: bytecount too large (more than 1073741822)
root [3]

Here is the modified version of the macro:

gen_Sparse3D.cxx (6.5 KB)

@pcanal May I request you to please help us to resolve this issue as @Wile_E_Coyote has also pointed out in the last post of this thread. It is very important in our data analysis to have a 3D-histogram to cleanly identify and assign gamma rays to a particular nucleus. Any help is greatly appreciated.

Thanks & regards,

Ajay