Draw and fit in a for loop

[newinthis]

Hello,

I have a TChain (called c) that have 100 TTrees. Each TTree store an image of 1280x1024 pixels, so it have 3 branches: x, y, pix_value.
To get the standard deviation of the value of each pixel, i am doing this:

TH1F* h0 = new TH1F("h0","h0",90,-10,10);
TH1F* hpn = new TH1F("hpn","hpn",100,6,10);
TF1 *noise = new TF1("noise","gaus",-10,10);

    for(Int_t x=0;x<1280;x++) {
      for(Int_t y=0;y<1024;y++) {
        sprintf(cuts,"x==%d && y==%d",x,y);
        c.Draw("pix_value>>h0",cuts,"goff"); 
        h0->Fit("noise","r");
        hpn->Fill(noise->GetParameter(2));
      }   
    }

This is taking infinity. Is there a correct way to do this?

Thank you!

[bellenot]

Try to use ACLiC (i.e. compiling your code), it might be faster…

[newinthis]

Thank you for the answer. I am compiling the code, and it also takes a lot of time.

Looks like the Draw function is slow, at least to run it 1millon of times.

[bellenot]

I suppose that both the draw and the fit are quite time consuming…

[Eddy_Offermann]

So you are accessing 100 1280x1024 images and would
like to get the jitter in in each of these pixels, resulting in

more than 1 million calls to create a histogram and then fit a
Gaussian line shape in order to get the spread …

For starters I would not fill in this very cpu-expensive way (define cut
a histogram just to get the (I,j) entry in the tree.

Why not fill an vector with the 100 (I,j) entries and just calculate
<x - > etc.

At some point in your analysis you created these trees. That could be the better entry point to calculate this noise.

-Eddy

[Triple_S]

Hi,
I have two proposals for you:

1. First project the tree into a TH3F and then do the loop. This should be already much faster, because projecting histograms is much less resource consuming than projecting trees. The code would be:

TH3F* hpixels = new TH3F("hpixels", "", 1280, -0.5, 1279.5, 1024, -0.5, 1023.5, 90, -10., 10.);
c.Project("hpixels", "x:y:pix_value");

TH1F* h0 = NULL; 
TH1F* hpn = new TH1F("hpn","hpn",100,6,10);
TF1 *noise = new TF1("noise","gaus",-10,10); 

for(Int_t x=0;x<1280;x++) {
    for(Int_t y=0;y<1024;y++) { 
        h0 = hpixels->ProjectionZ("h0", x+1, x+1, y+1, y+1);
        h0->Fit("noise","r"); 
        hpn->Fill(noise->GetParameter(2));
        delete h0;
    } 
}

2. Since you said you only want to get the standard deviation of each pixel, why not calculate it manually with its definition σ = sqrt(sum for i=1 to N (x - μ)^2/(N-1)) of course you need the mean values, but these can be calculated already by two 2D histograms, see the following code:

TH1F* hpn = new TH1F("hpn","hpn",100,6,10);

TH2F* hMean = new TH2F("hMean", "", 1280, -0.5, 1279.5, 1024, -0.5, 1023.5);
TH2F* hEntries = new TH2F("hEntries", "", 1280, -0.5, 1279.5, 1024, -0.5, 1023.5);
c.Project("hMean", "x:y", "pix_value");
c.Project("hEntries", "x:y");
hMean->Divide(hEntries);

Int_t x, y;
Float_t pix_value;
c.SetBranchAddress("x", &x);
c.SetBranchAddress("y", &y);
c.SetBranchAddress("pix_value", &pix_value);

TH2F* hSigmaSum = new TH2F("hSigmaSum", "", 1280, -0.5, 1279.5, 1024, -0.5, 1023.5);
for (Long64_t i = 0; i < c.GetEntries(); i++) {
    c.GetEntry(i);
    hSigmaSum->Fill(x, y, TMath::Power(pix_value - hMean->GetBinContent(x+1, y+1), 2.));
}

for (x = 0; x < 1280; x++)
    for(y = 0; y < 1024; y++)
        hpn->Fill(TMath::Sqrt(hSigmaSum->GetBinContent(x+1, y+1) / (hEntries->GetBinContent(x+1, y+1) - 1)));

Both codes should produce the same results but be much faster than your version. I have not tested it, but I think the second one should be the fastest.

[Eddy_Offermann]

I just want to add to @Triple_S suggestions that in general using histograms to get the variances is overkill.

Assuming data array x:

double x_sum = 0.;
double xSq_sum = 0.;
int n =0;
std::vector x(..);
for (..) {
   x_sum += x[i];
   xSq_sum += x[i]*x[i];
   n++;
}
double x_var = (xSq_sum-std::pow(x_sum,2))/(n-1);