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

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))

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.

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.

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

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

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)

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
}

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

I think you can easily reduce the amount of the required RAM by (up to) a factor 6, if you explicitly sort your gammas:

#include <algorithm>
// ...
      for(int j = 0; j < (No_Clover - 2); j++) {
        for(int k = (j + 1); k < (No_Clover - 1); k++) { // (k > j)
          for(int l = (k + 1); l < No_Clover; l++) { // (l > k)
            // ...
                        // Double_t x[3] = {(Double_T_t)Eclab[j], (Double_T_t)Eclab[k], (Double_T_t)Eclab[l]};
                        // Double_t x[3] = {(Eclab[j] + 1.e-6), (Eclab[k] + 1.e-6), (Eclab[l] + 1.e-6)};
                        Double_t x[3] = {(Eclab[j] + 0.5), (Eclab[k] + 0.5), (Eclab[l] + 0.5)};
                        std::sort(x, x + 3); // sort in an increasing order
                        if (his3D->GetNbins() <= MaxNbins) his3D->Fill(x);
                        // ...

When you do not sort your gammas, each triple of energies appears 6 times in your 3D histogram so, there are 6 “peaks” in different areas of the 3D histogram (where 6 is the number of permutations, of course).

If you sort them (in an increasing order), you will get 1 “peak” only (always energy_gamma_1 <= energy_gamma_2 <= energy_gamma_3).

This makes your 3D histogram “totally unsymmetric” but, I think you can take it into account when projecting / analysing it later (well, I’m afraid that this idea requires some detailed thoughtful examination so, you should discuss it with your colleagues).

1 Like

I am open to make thorough checks using the procedure that you have mentioned. However, I don’t to how to analyze such “unsymmetric” cube. Can you please tell me the possible analysis method/s?

If I have to make such a “totally unsymmetric” cube, then probably, TH3F class will also work. Then, I can make use of TSpectrum3 for subtracting 3D-background from the 3D-histogram, which I guess not possible in THnSparse.

BTW, I am able to successfully generate 3D THnSparse with lower dimensions (1000) with full data set.

Is this same as scaling or using weightage while filling the histogram?
i.e. his3D->Fill(x,0.1666666);

Maybe this works:

#include "TH1.h"
#include "TH2.h"
#include "THnSparse.h"
#include "TAxis.h"
#include "TString.h"

// Usage example: TH1D *h1 = Project1D(his3D, 449.5, 452.5, 668.5, 671.5);
// Note: in principle, you should set low1 <= up1 and low2 <= up2
TH1D *Project1D(THnSparse *h, // the ("totally unsymmetric") "cubic" histogram
                Double_t low1, Double_t up1, // the first "window"
                Double_t low2, Double_t up2, // the second "window"
                Option_t *option = "") { // THnSparse::Projection option
  if (!h) return 0; // just a precaution
  TAxis *a0 = h->GetAxis(0);
  TAxis *a1 = h->GetAxis(1);
  TAxis *a2 = h->GetAxis(2);
  a0->SetRange(a0->FindFixBin( low1 ), a0->FindFixBin( up1 ));
  a1->SetRange(a1->FindFixBin( low2 ), a1->FindFixBin( up2 ));
  a2->SetRange(-1, -1);
  TH1D *h1 = h->Projection(2, option);
  h1->SetName(TString::Format("%s_Project1D", h->GetName()));
  a0->SetRange(a0->FindFixBin( low2 ), a0->FindFixBin( up2 ));
  a1->SetRange(a1->FindFixBin( low1 ), a1->FindFixBin( up1 ));
  TH1D *hi = h->Projection(2, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(a0->FindFixBin( low1 ), a0->FindFixBin( up1 ));
  a1->SetRange(-1, -1);
  a2->SetRange(a2->FindFixBin( low2 ), a2->FindFixBin( up2 ));
  hi = h->Projection(1, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(a0->FindFixBin( low2 ), a0->FindFixBin( up2 ));
  a2->SetRange(a2->FindFixBin( low1 ), a2->FindFixBin( up1 ));
  hi = h->Projection(1, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(a1->FindFixBin( low1 ), a1->FindFixBin( up1 ));
  a2->SetRange(a2->FindFixBin( low2 ), a2->FindFixBin( up2 ));
  hi = h->Projection(0, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a1->SetRange(a1->FindFixBin( low2 ), a1->FindFixBin( up2 ));
  a2->SetRange(a2->FindFixBin( low1 ), a2->FindFixBin( up1 ));
  hi = h->Projection(0, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);
  return h1;
}

// Usage example: TH2D *h2 = Project2D(his3D, 449.5, 452.5);
// Note: in principle, you should set low <= up
TH2D *Project2D(THnSparse *h, // the ("totally unsymmetric") "cubic" histogram
                Double_t low, Double_t up, // the "window"
                Option_t *option = "") { // THnSparse::Projection option
  if (!h) return 0; // just a precaution
  TAxis *a0 = h->GetAxis(0);
  TAxis *a1 = h->GetAxis(1);
  TAxis *a2 = h->GetAxis(2);
  a0->SetRange(a0->FindFixBin( low ), a0->FindFixBin( up ));
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);
  TH2D *h2 = h->Projection(1, 2, option);
  h2->SetName(TString::Format("%s_Project2D", h->GetName()));
  TH2D *hi = h->Projection(2, 1, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(a1->FindFixBin( low ), a1->FindFixBin( up ));
  a2->SetRange(-1, -1);
  hi = h->Projection(0, 2, option);
  h2->Add(hi);
  delete hi; // no longer needed
  hi = h->Projection(2, 0, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(a2->FindFixBin( low ), a2->FindFixBin( up ));
  hi = h->Projection(0, 1, option);
  h2->Add(hi);
  delete hi; // no longer needed
  hi = h->Projection(1, 0, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);
  return h2;
}
1 Like

Please give me some time to test this. Will get back to you with the outcome after testing.

In the meantime, can you please answer my last question?

You have mentioned in the post 61 of this thread that we can reduce the RAM usage by a factor of 6 by creating an “asymmetric” 3D-histogram. Is it the same as using

his3D->Fill(x,0.1666666); // 1/6 = 0.16666666

instead of

his3D->Fill(x);

Just try it.

I tried it, and it failed to work :cry:

Now, I am testing the last method. Will get back to you with results.

BTW, do you know whether @pcanal is working on/has time to resolve the issue?

I have tested the “unsymmetric” cubic histogram and the corresponding projection method that you have suggested. The preliminary observations indicate that everything is working perfectly as expected. I will do more thorough check in the weekend.

As I had mentioned earlier, I have to create several projection1D using different combination of gamma-ray pairs. For this, I am using the macro “Sparse_project.C” (attached). I use the following sequence to generate two 1D-Projections:

root [0] TFile *f = TFile::Open("sparse3D_asym_all_4k.root");
root [1] .ls
TFile**		sparse3D_asym_all_4k.root	
 TFile*		sparse3D_asym_all_4k.root	
  KEY: THnSparseT<TArrayF>	sparse3D;1	sparse3D
root [2] .L Sparse_project.C 
root [3] Project1D(sparse3D,449.5,453.5,668.5,672.5);
root [4] .ls
TFile**		sparse3D_asym_all_4k.root	
 TFile*		sparse3D_asym_all_4k.root	
  OBJ: TH1D	sparse3D_Project1D	sparse3D projection EclabZ : 0 at: 0xec6d440
  KEY: THnSparseT<TArrayF>	sparse3D;1	sparse3D
root [5] sparse3D_Project1D->Draw();
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
root [6] sparse3D_Project1D->SetName("his_451_670");
root [7] his_451_670->Draw();
root [8] Project1D(sparse3D,289.5,293.5,698.5,702.5);
root [9] .ls
TFile**		sparse3D_asym_all_4k.root	
 TFile*		sparse3D_asym_all_4k.root	
  OBJ: TH1D	his_451_670	sparse3D projection EclabZ : 0 at: 0xec6d440
  OBJ: TH1D	sparse3D_Project1D	sparse3D projection EclabZ : 0 at: 0xf129740
  KEY: THnSparseT<TArrayF>	sparse3D;1	sparse3D
root [10] sparse3D_Project1D->Draw();
root [11] sparse3D_Project1D->SetName("his_291_700");
root [12] his_291_700->Draw();
root [13] his_451_670->Draw();
root [14]

However, there is a problem. Both the spectra are identical. As you can see, I am using entirely different sets of gamma rays to generate two 1D-histograms. Am I doing anything wrong?

Also, I feel that it would be ideal to mention only gamma-ray energies while projecting, along with a variable “width parameter” which will internally set the “window” for the gamma rays.

Sparse_project.C (3.4 KB)

gen_Sparse3D.cxx (6.9 KB)

Try (first you need to root [0] .L Sparse_project.C, of course):

{
  TFile *f = TFile::Open("sparse3D_asym_all_4k.root");
  THnSparse *h; f->GetObject("sparse3D", h);
  Double_t g1 = 451.5;
  Double_t g2 = 670.5;
  Double_t w = 2.0;
  TH1D *h1 = Project1D(h, g1 - w, g1 + w, g2 - w, g2 + w);
  h1->SetName(TString::Format("his_%d_%d", Int_t(g1), Int_t(g2)));
  h1->Print();
  g1 = 291.5;
  g2 = 700.5;
  TH1D *h2 = Project1D(h, g1 - w, g1 + w, g2 - w, g2 + w);
  h2->SetName(TString::Format("his_%d_%d", Int_t(g1), Int_t(g2)));
  h2->Print();
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  h1->Draw();
  c->cd(2);
  h2->Draw();
  c->cd(0);
  // do we get exactly the same output as above?
  h1->Print();
  h2->Print();
}

It worked, thank you. The two 1D-histograms are no longer identical.

Please check the attached modified version of the macro which I used for 1D- and 2D- projections of THnSparse.

We know that we have TSpectrum2 and TSpectrum3 classes available for background generation of TH2 and TH3 histograms, respectively.

Are there any equivalent “background generation” methods/procedures available for THnSparse Class?

Sparse_project.C (3.7 KB)

You don’t need to care if “eg_low <= eg_hi”.
Use “SetName” once per histogram (just modify the first occurrence).

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