Filling a 3 dimensional array of histograms

Hello, I am handling a large amount of data with a 3 dimensional array of histograms with indices A, P and M. A goes from 0 to 2, P goes from 0 to 4 and M from 0 to 11, and for some reason, I can store things in the coordinates of the form [0][P][M] but not in the coordinates of the form [1][P][M] and [2][P][M].

Here is the code:

First, some functions I use:

template<typename T, int N>
int getSize(T (&arr1)[N]) //Passing the array by reference 
{
   return N; //Correctly prints the size too [cool trick ;-)]
}

int getRange(double v, double arr[], int n)
{
  for (int i = 0; i < n; ++i)
  {
    if (v>=arr[i])
    {
      return i;
    }
  }
}

int getMultRange(Event & event, double mult_range[], int n)
{
  double mult=0; 
  for (int i = 0; i <= event.size(); ++i) if (event[i].isCharged() && event[i].isFinal())
  {
    ++mult;
  }
  return getRange( mult, mult_range, n);
}

The purpose for getRange and getMultRange is because in my code I often want to find out in what point in the distribution a certain value falls, so I have this vectors with the distribution ranges from higher to lower.

Now the relevant code

    double pT_range[]         = {8, 4, 2, 1, 0.5}; //multiplicity ranges
    double Aj_range[]         = {0.22, 0, 0};
    double mult_range[]       = {234, 204, 188, 168, 152, 140, 128, 116, 104, 92, 76, 0}; //multiplicity ranges
    TH2D*  all_histos[getSize(Aj_range)][getSize(pT_range)][getSize(mult_range)];
    const  Char_t* forHist_Aj_range[] = {"AjInc", "Aj>0.22", "Aj<0.22"};
  for (int a = 0; a < getSize(Aj_range); ++a)
  {
    for (int p = 0; p < getSize(pT_range); ++p)
    {
      for (int m = 0; m < getSize(mult_range); ++m)
      {
        all_histos[a][p][m] = new TH2D( Form("hist_mult>=%f_Aj%s_p>=%f",mult_range[m],forHist_Aj_range[a],pT_range[p]), "The momentum of particles in each Delta Ring", 9, 0, 1.8, 1000, -20, 20);
      }
    }
  }
//I've also tried with:
//TH2D*  all_histos[3][5][12];
//... in between we get the leading and subleading jets of the event    
    Aj = (leading.perp()-subleading.perp())/(leading.perp()+subleading.perp()); // we get the centrality
    a = getRange( Aj, Aj_range, getSize(Aj_range)) + 1; //we determine what range of centrality this is
    m = getMultRange( event, mult_range, getSize(mult_range)); //we determine what range of multiplicity this is
    for (int n = 0; n <= event.size(); ++n) if (event[n].isCharged() && event[n].isFinal())
    {
      particle = event[n]; //we declare a particle
      delta = sqrt(pow(particle.eta()-axis.eta(),2)+pow(particle.phi()-axis.phi(),2)); //we calculate the delta that we want to store in the histogram

      if ( cutContidion(particle, delta) ) continue; //cutCondition is a custom function 
      pTll = -particle.perp()*cos( particle.phi()-axis.phi() ); //we calculate the "projection" that we want to store in the histogram

      p = getRange( particle.perp(), pT_range, getSize(pT_range)); //we determine what range of transverse momentum this is

      all_histos[0][p][m]->Fill(delta,pTll); //to fill Aj Inclusive
cout<<"This will be printed"<<endl;
      all_histos[a][p][m]->Fill(delta,pTll); //to fill either Aj>0.22 or Aj<0.22 
cout<<"This will NOT be printed, the code fails in the previous line"<<endl;
    }

The code compiles, but the error I get when trying to run this code is:

*** Break *** segmentation violation
===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0 0x0000003703aac7be in waitpid () from /lib64/libc.so.6
#1 0x0000003703a3e5c9 in do_system () from /lib64/libc.so.6
#2 0x00007f35f722b708 in TUnixSystem::StackTrace (this=0x103b880) at /storage/icn/aortizve/software/64bits/alicesw/root/alice_v5-34-30/src/core/unix/src/TUnixSystem.cxx:2419
#3 0x00007f35f722ab93 in TUnixSystem::DispatchSignals (this=0x103b880, sig=kSigSegmentationViolation) at /storage/icn/aortizve/software/64bits/alicesw/root/alice_v5-34-30/src/core/unix/src/TUnixSystem.cxx:1294
#4
#5 0x000000000046e96c in main ()
===========================================================
The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
Sign in to GitHub · GitHub. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5 0x000000000046e96c in main ()
===========================================================

I am at the point where I can’t imagine what the problem may be, I’ve checked that the indices are within bounds, that the indices are “int”, if I never try to store data in those histograms the code runs and if I save all_histos the histograms are indeed there, although empty.

Thanks a lot.

Edit: I solved it?

My code works now, but while fixing it I found that it shouldn’t have worked in the first place. The mistake is almost embarrassing: pT was being smaller than 0.5, the low boundary of the last considered interval. So I added to cutCondition an extra condition to remove those events but then I realized that it WAS filling the histograms with index [0][P][M], which is weird, if anything the problem was with the index P, but it didn’t “manifest” until the index A=1,2, and those histograms had data I could see it…

But whatever

1 Like

You don’t seem to create these histograms anywhere.

Forgot to mention that I use this, I edited my post to include them:

  const  Char_t* forHist_Aj_range[] = {"AjInc", "Aj>0.22", "Aj<0.22"};
  for (int a = 0; a < getSize(Aj_range); ++a)
  {
    for (int p = 0; p < getSize(pT_range); ++p)
    {
      for (int m = 0; m < getSize(mult_range); ++m)
      {
        all_histos[a][p][m] = new TH2D( Form("hist_mult>=%f_Aj%s_p>=%f",mult_range[m],forHist_Aj_range[a],pT_range[p]), "The momentum of particles in each Delta Ring", 9, 0, 1.8, 1000, -20, 20);
      }
    }
  }

The easiest thing could be to add (right before the line which fails):

std::cout << a << " " << p << " " << m << std::endl;

Try also:

and:

1 Like

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