# Draw and fit in a for loop

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!

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

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.

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

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

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;
}
}
``````
1. 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;

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.

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

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