Saturday, August 09, 2014

data analysis with ROOT - redux

In my previous post about data analysis with ROOT, I had used the data from a single gamma. Extending it to multiple gammas - 1500 - and calculating the energy weighted position, found a bug in my code.

Apparently the Fill() function does not fill if the indices passed to it are 0. That is,
h2->Fill(i,j,energy2D[i][j])
will not work for i=0 or j=0. Since our energy2D C array has a range 0-7 for i & j, this Fill() had to be re-written as
h2->Fill(i,j,energy2D[i-1][j-1])
with the i & j range 1-8.

I also changed my definition of the histogram to
TH2I *h2 = new TH2I("h2","energy distribution",8,1,9, 8,1,9);
instead of
TH2I *h2 = new TH2I("h2","energy distribution",8,0,8, 8,0,8);

With the resulting root file, energy-weighted position analysis could be done quite quickly, almost instantaneously, with

// loading the data in root file
    TFile f("multifile.root");
    f.ls();
    TH2I* h = (TH2I*)f.Get("h2"); 

for (int i=1;i<9;i++)
    {     
    for (int j=1;j<9;j++)
      {
      bincontent=h->GetBinContent(i,j);
      ewposx += posx[i-1][j-1] * bincontent ;
      ewposy += posy[i-1][j-1] * bincontent ;
      totalcount +=  bincontent ;      
      }
    }
    ewposx /= totalcount;
    ewposy /= totalcount;
    printf("ewposx & y   = %f, %f\n", ewposx  , ewposy  );

Changing the position of the gamma was fairly straightforward in LXePrimaryGeneratorAction.cc

The calculation of energy weighted position in the LXe example seems to have broken somewhere in development - people have reported non-zero positions in 2010 and so on, but when it is run now, it returns only zero, since the GetPos() for the pmt hit collection class returns zero. This seems to have causes in the class definitions and so on - LXePMTHit.cc doesn't mention fPos private variable in the over-ride function for = operator and so on. Instead of trying to fix that, I just re-entered the position data in root, as

for (int i=0;i<8;i++)
    {
    y = -22.75;
    for (int j=0;j<8;j++)
      {
      posx[i][j]= x ;
      posy[i][j]= y ;

      //printf("pos2D[%d,%d]=%f,%f\n",i,j,posx[i][j] , posy[i][j] ); 
      y+=6.5;
      }
    x+=6.5;
    }



No comments:

Post a Comment