Multiple template fitting a pulse

[benkrikler]

So I'm moving on with fitting multiple templates to a TPI in the case of pile-up. It's not as straight forward as it first sounds though.

The main question I have right now is about forming the chi-2: When I come to fit the pile-up with 2 templates i will only fit the window where all templates overlap (plus some safety ). Otherwise, I would need to extrapolate the undershoot of the first template somehow to know the baseline for the second pulse. When I then loop over this window for a given set of fit parameters (time offsets, amplitude factors and pedestal) for each bin I:

1. calculate the residual for the fitted templates (the sum of all template*scaling - actual TPI sample value).
2. Divide the square of this by the square of the templates' error for that sample.

So my question is, what is now the error for each sample in the template? My immediate feeling is to just add in quadrature the corresponding samples' errors from all templates. Since the error bars on the template are the RMS of all samples that contributed to that bin this felt ok, but I wanted tose what others thought.

[benkrikler]

I forgot the scaling factors need to go into the errors, so the square of the total error for a given sample would be the sum over all templates of the product of the scaling factor and the square of the template's error for that sample. (Maths is hard in words!)

Total error ^2= scale factor1 * ( template error 1) ^2 + scale factor2 * ( template error 2) ^2

[benkrikler]

Actually, that might explain the crazy big chi-2 I showed for the fits today, since I dont think the chi-2 calculation incorporates the scale factor for the weights. That would give a correction of N* scale_factor where N is the number of samples in the fit since it would appear in the dominator for each term in the sum which gives the chi-2. With roughly 200 samples and an amplitude of around 800 as in the pulses I showed earlier, that would give a factor of 1.6e5 which would look almost spot on!

[litchfld]

I forgot the scaling factors need to go into the errors, so the square of the total error for a given sample would be the sum over all templates of the product of the scaling factor and the square of the template's error for that sample. (Maths is hard in words!)

Total error ^2= scale factor1 * ( template error 1) ^2 + scale factor2 * ( template error 2) ^2

You mean:

Yes I think that is sensible - for the template, see below. The errors when forming templates is a bit fiddly because we are both iterating and trying to estimate the error from the random noise overlay, but for because of that it should be sensible to apply them in the obvious way. (BTW I assume the errors are scaled down in the same way when the templates are normalised? Otherwise it doesen't make sense.)

I'm not 100% sure whether there should be an additional contribution from the pulse. I think not because it's dominated by random noise, and that's the error we try to account for when making the template. But there might be a small additional contribution that we aren't accounting for from uncorrelated bin-wise fluctuations? Try it and see.

[benkrikler]

I don't seem to be able to see your formula their Phill, could you upload it again?

The errors when forming templates is a bit fiddly because we are both iterating and trying to estimate the error from the random noise overlay

I don't really understand what you're saying here. You mean because we're scaling each pulse to fit into the template it's not clear what error to add in with each iteration?

BTW I assume the errors are scaled down in the same way when the templates are normalised? Otherwise it doesen't make sense.

Actually, that's a very good point, I don't think I'd checked that properly. This is the heart of the normalisation code as it is:

    // Subtract off the pedesal
for (int iBin = 1; iBin <= fTemplatePulse->GetNbinsX(); ++iBin) {
  double old_value = fTemplatePulse->GetBinContent(iBin);
  double new_value = old_value - template_pedestal;

  fTemplatePulse->SetBinContent(iBin, new_value);
}

double norm = std::fabs(fTemplatePulse->GetMaximum());  // for positive polarity templates
fTemplatePulse->Scale(1.0/norm);

So all the work to scale the errors is buried in TH1::Scale. This is a bit of ROOT I'm not very familiar with though. I believe we need to declare the template an average type histogram somehow and then call SumW2 before we call scale. Have I got that right? If we're not sure, i can just do it by hand in the loop above.

[AndrewEdmonds11]

From the ROOT manual, it looks like you're right mostly right. I'm not sure what an "average" histogram is, but here it just seems to say that we should call Sumw2 first.

void Scale(Double_t c1 = 1, Option_t* option = "")
  -*-*-*Multiply this histogram by a constant c1*-*-*-*-*-*-*-*-*

   this = c1*this

 Note that both contents and errors(if any) are scaled.
 This function uses the services of TH1::Add

 IMPORTANT NOTE: If you intend to use the errors of this histogram later
 you should call Sumw2 before making this operation.
 This is particularly important if you fit the histogram after TH1::Scale

 One can scale an histogram such that the bins integral is equal to
 the normalization parameter via TH1::Scale(Double_t norm), where norm
 is the desired normalization divided by the integral of the histogram.

 If option contains "width" the bin contents and errors are divided
 by the bin width.

[litchfld]

The errors when forming templates is a bit fiddly because we are both iterating and trying to estimate the error from the random noise overlay

I don't really understand what you're saying here. You mean because we're scaling each pulse to fit into the template it's not clear what error to add in with each iteration?

Not exactly. Because you don't have any initial estimate of the random errors on each individual pulse ( sqrt(n) is clearly not right) you need to iterate in such a way that the errors are also estimated, and also converge. I think the iterated error that comes with the template will cover everything, including random variations on any pulse you subsequently 'measure' with the template. But i'm not 100% sure off the top of my head. My point was just that the actual measurement using the finished template should be much more straightforward, so this shouldn't need over thinking...

[litchfld]

This is the heart of the normalisation code as it is:

for (int iBin = 1; iBin <= fTemplatePulse->GetNbinsX(); ++iBin) {
  //...
  fTemplatePulse->SetBinContent(iBin, new_value);
}
//...
fTemplatePulse->Scale(1.0/norm);

If that everything then i think we are fine. The function SetBinContent() will not modify the errors so they should remain from the construction, and Scale will apply the same factor as to the bin content.

[thnam]

Hi Phil, Please have a look at the EdE tree if you can: https://www.dropbox.com/s/gtf4t0qgj3b0pdq/si16p_edetree.root?dl=0 The file is nearly 400 MB …

The energy of silicons; timing of muSc and thick silicons are stored. To make dE/E plot for 1 quadrant: ede_tree->Draw(“e_sir11:e_sir11+e_sir2”)

Cheers, Nam=

[litchfld]

Got it thanks... why did the github link take me here??