Measuring Darkcount rate

dneise commented 6 years ago

In PR #175 @calispac wrote this:

Maybe to have better dark count rate measurement one could use Poisson count method (i.e. number of zero events). In this algorithm (spe.py) the zeros are not counted since we explicitly ask for finding pulses (above a certain thershold ~ 0.5 p.e.). What was told to me by @AndriiNagai is that we could just make a histogram of (np.max(adc_samples, axis=-1) - np.min(adc_count, axis=-1)).

This will not be so good for gain measurement since the 0. p.e peak is shifted to higher ADC count but would allow us to retrieve dark count rate. Also I doubt that the fit will work with the 0 p.e. peak.

To create such a histogram as the code is, we need to implement a find_pulse_x() that fills pulse_mask with np.arange(0, adc_samples.shape[-1]) == np.argmax(adc_samples, axis=-1) and have a function that fills the baseline Field with np.min(adc_samples, axis=-1)

And I thought this discussion is worth its own issue.

dneise commented 6 years ago

@AndriiNagai is completely right that using the so called "Poisson method" is a better way to estimate the dark count rate.

I however have not yet heard of this particular implementation of this method. But I guess everything is worth a try.

AndriiNagai commented 6 years ago

Hi Dominik,

On 8 Apr 2018, at 18:30, Dominik Neise notifications@github.com wrote:

@AndriiNagai https://github.com/AndriiNagai is completely right that using the so called "Poisson method" is a better way to estimate the dark count rate.

I however have not yet heard of this particular implementation of this method. But I guess everything is worth a try.

I have used this method during my PhD. Please, find attached the comparison between DCR calculated from Poisson method and from measurements with Oscilloscope, where threshold was set to 0.5 p.e amplitude and time difference between two trigged events was used to calculate DCR. You can see good agreement between both methods at T = 258.15K, while at higher T measurements with Oscilloscope saturated due to Oscilloscope dead time

Bets regards, Andrii

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379563101, or mute the thread https://github.com/notifications/unsubscribe-auth/AYRh4FR2z2-Y8Bhir19iz289qKyEEwlLks5tmjsVgaJpZM4TLmcT.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/cta-sst-1m/digicampipe","title":"cta-sst-1m/digicampipe","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/cta-sst-1m/digicampipe"}},"updates":{"snippets":[{"icon":"PERSON","message":"@dneise in #176: @AndriiNagai is completely right that using the so called \"Poisson method\" is a better way to estimate the dark count rate. \r\n\r\nI however have not yet heard of this particular implementation of this method. But I guess everything is worth a try."}],"action":{"name":"View Issue","url":"https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379563101"}}}

AndriiNagai commented 6 years ago

This Poisson method is widely used during absolute PDE measurements. To correct the PDE to DCR, please see attached (Section 2.2).

Bets regards, Andrii

On 9 Apr 2018, at 09:47, Andrii Nagai Andrii.Nagai@unige.ch wrote:

Hi Dominik,

On 8 Apr 2018, at 18:30, Dominik Neise <notifications@github.com mailto:notifications@github.com> wrote:

@AndriiNagai https://github.com/AndriiNagai is completely right that using the so called "Poisson method" is a better way to estimate the dark count rate.

I however have not yet heard of this particular implementation of this method. But I guess everything is worth a try.

I have used this method during my PhD. Please, find attached the comparison between DCR calculated from Poisson method and from measurements with Oscilloscope, where threshold was set to 0.5 p.e amplitude and time difference between two trigged events was used to calculate DCR. You can see good agreement between both methods at T = 258.15K, while at higher T measurements with Oscilloscope saturated due to Oscilloscope dead time

Bets regards, Andrii <Screen Shot 2018-04-09 at 09.40.28.png>

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379563101, or mute the thread https://github.com/notifications/unsubscribe-auth/AYRh4FR2z2-Y8Bhir19iz289qKyEEwlLks5tmjsVgaJpZM4TLmcT.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/cta-sst-1m/digicampipe","title":"cta-sst-1m/digicampipe","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/cta-sst-1m/digicampipe https://github.com/cta-sst-1m/digicampipe"}},"updates":{"snippets":[{"icon":"PERSON","message":"@dneise in #176: @AndriiNagai is completely right that using the so called \"Poisson method\" is a better way to estimate the dark count rate. \r\n\r\nI however have not yet heard of this particular implementation of this method. But I guess everything is worth a try."}],"action":{"name":"View Issue","url":"https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379563101 https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379563101"}}}

dneise commented 6 years ago

@AndriiNagai wrote:

This Poisson method is widely used during absolute PDE measurements

Yes .. I absolutely agree. The Poisson method is widely used. But this particular implementation I have not heard of. By particular implementation, I mean using the difference between min and max of a sampled timeline. I usually see the Poisson method being applied to a histogram which is the result of some kind charge extraction (usually integrating in a small time window around some kind of trigger, random trigger is also possible here).

That's all I'm saying. I did not want to criticize anything here. Neither the method, which is indeed well known, I think I saw it first 10 years ago used by Dieter Renker (the author of this very nice overview article ..Advances in solid state photon detectors, but I am not sure where this method was used the first time ... I assume it is much older.

However, Dieter always integrated a bit, and did not use max-min.

The attachment unfortunately did not make it from the mail into this thread here. Maybe you could give us a link?

AndriiNagai commented 6 years ago

On 9 Apr 2018, at 10:30, Dominik Neise notifications@github.com wrote:

@AndriiNagai https://github.com/AndriiNagai wrote:

This Poisson method is widely used during absolute PDE measurements

Yes .. I absolutely agree. The Poisson method is widely used. But this particular implementation I have not heard of. By particular implementation, I mean using the difference between min and max a sampled timeline. I usually see the Poisson method being applied to a histogram which is the result of some kind charge extraction (usually integrating in a small time window around some kind of trigger, random trigger is also possible here).

That's all I'm saying. I did not want to criticize anything here. Neither the method, which is indeed well known, I think I saw it first 10 years ago used by Dieter Renker (the authot of this very nice overview article ..Advances in solid state photon detectors http://iopscience.iop.org/article/10.1088/1748-0221/4/04/P04004/pdf, but I am not sure where this method was used the first time ... I assume it is much older.

However, Dieter always integrated a bit, and did not use max-min.

For commercial SiPM devices, at a given T and Vbias, charge is proportional to amplitude. So, for me, between “integral" and "max-min" is only the difference that "max-min" (I would better write "max-baseline") is easier to calculate, because integral is more sensitive to baseline calculation. The attachment unfortunately did not make it from the mail into this thread here. Maybe you could give us a link?

Sure: https://www.sciencedirect.com/science/article/pii/S0168900210008156 https://www.sciencedirect.com/science/article/pii/S0168900210008156

May be, we can discuss this during the meeting?

Best regards, Andrii

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379675508, or mute the thread https://github.com/notifications/unsubscribe-auth/AYRh4CUBtY1KoYpuL17alHF5ZLbmOHpzks5tmxwtgaJpZM4TLmcT.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/cta-sst-1m/digicampipe","title":"cta-sst-1m/digicampipe","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/cta-sst-1m/digicampipe"}},"updates":{"snippets":[{"icon":"PERSON","message":"@dneise in #176: @AndriiNagai wrote:\r\n\u003e This Poisson method is widely used during absolute PDE measurements\r\n\r\nYes .. I absolutely agree. The Poisson method is widely used. But this particular implementation I have not heard of. By particular implementation, I mean using the difference between min and max a sampled timeline. I usually see the Poisson method being applied to a histogram which is the result of some kind charge extraction (usually integrating in a small time window around some kind of trigger, random trigger is also possible here).\r\n\r\nThat's all I'm saying. I did not want to criticize anything here. Neither the method, which is indeed well known, I think I saw it first 10 years ago used by Dieter Renker (the authot of this very nice overview article ..Advances in solid state photon detectors, but I am not sure where this method was used the first time ... I assume it is much older.\r\n\r\nHowever, Dieter always integrated a bit, and did not use max-min.\r\n\r\nThe attachment unfortunately did not make it from the mail into this thread here. Maybe you could give us a link?"}],"action":{"name":"View Issue","url":"https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-379675508"}}}

dneise commented 6 years ago

@AndriiNagai this "max-min" method seems indeed to work nicely.

I tested it on toy monte-carlos (3 MHz dark count rate) against the method I had in mind, which is integrating at random places.

(NB: I seem to have a systematic error of some kind still .. I am over-estimating with both methods) But still both methods clearly estimate something around 3MHz and the max-min method is faster in execution.

I also applied this method to real data ("SST1M01_20171124.001.fits.fz") once:

Here also both methods seem to estimate something similar ...

calispac commented 6 years ago

Awesome ! Could you tell me what is random charge method ? Also If we can merge what you have done in #182 that would be nice.

dneise commented 6 years ago

Could you tell me what is random charge method ?

Sure sorry.

Imagine this lab setup:

you have pulsed light source with a very low expectation value for the number of photons, e.g. 0.1 photon per shot.
you have your PMT/SiPM, a pre-amp and a normal oscilloscope.
you trigger the scope with the trigger from the light source.
at the trigger time you integrate the signal on the scope for ~5ns or so (whatever the typical pulse width is)

Then you get a similar spectrum as we get with our SPE-analysis, but since your light source is so weak, you will very often get a trigger, but no pulse at all, resulting in a well visible 0-p.e. peak in the spectrum.

From this 0-p.e. peak in the spectrum one can nicely estimate the expectation value of the light source.

Now with dark counts, we do not have a trigger, but we can simply think of dark counts as a weak light source. In any given integration window (of e.g. 10ns width) the probability of having a dark count is very low (at 3MHz it is 3%). So we can simply create a random trigger, and make the same measurement as the one we did in the lab with a triggered light source.

This is basically the "Poisson method" as far as I understood it.

Now the analysis goes like this:

count the number of entries in the 0-p.e. peak: N0
count the total number of entries. N
frequentist probability for "no pulse in time window" P0 = N0 / N
Poisson probability for "exactly zero events in window" = exp(- mu)
where mu is the expectation value: mu = time_window * average_rate

So one can calculate the rate as:

average_rate = - ln(N0/N) / time_window

calispac commented 6 years ago

Ok so this means that we are underestimating the N0 if the dark count rate is overestimated ?

This seems normal if the 0 p.e. peak and the 1 p.e. peak are not well separated. See an example of this (max- baseline) histogram pixel_8

(By the way I used most probable value for baseline not min)

dneise commented 6 years ago

Yes, underestimating N0 is one possibility.

How was this spectrum made?

AndriiNagai commented 6 years ago

Nice work Dominik!

From my experience, the tricky point is to separate 0p.e. from 1p.e. peak. Typically, the histogram contains entries between 0 p.e. and 1 p.e. how did the algorithm decide if entry with 0.5 p.e amplitude go to N0 or not? Could you show the example of max - min distribution? Is it possible to use something like: “max - baseline”, where baseline is an average over given time window before t_max - SiPM_risetime?

Another, very stupid question: Are you sure that you generated exactly 3MHz of DCR? Not 3.1 or 2.9? Could your waveform contain only part of SiPM pulses (if pulse was generated in the beginning or in the end of waveform)?

Best regards, Andrii

On 10 Apr 2018, at 17:09, Dominik Neise notifications@github.com wrote:

Yes, underestimating N0 is one possibility or overestimating N.

How was this spectrum made?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-380136255, or mute the thread https://github.com/notifications/unsubscribe-auth/AYRh4Chn3bj-9o_W3R4ucc4K5T1FMGYEks5tnMstgaJpZM4TLmcT.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/cta-sst-1m/digicampipe","title":"cta-sst-1m/digicampipe","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/cta-sst-1m/digicampipe"}},"updates":{"snippets":[{"icon":"PERSON","message":"@dneise in #176: Yes, underestimating N0 is one possibility or overestimating N.\r\n\r\nHow was this spectrum made? "}],"action":{"name":"View Issue","url":"https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-380136255"}}}

dneise commented 6 years ago

Unortunately I have no time to further go into details with this study at the moment. But all the details can be found here.

https://github.com/dneise/study_cta_sst1m_darkcount_rate

I might find time to add a readme.

calispac commented 6 years ago

average_rate = - ln(N0/N) / time_window

@dneise I just realized that the time_window is not 50 ns (e.g. time of readout window) but rather an effective window on which pulses are searched for. In this function : https://github.com/cta-sst-1m/digicampipe/blob/4ab7ae2c22387a0d9c0ec5214d03d262e1703631/digicampipe/scripts/spe.py#L215-L237

Pulses are searched over : T = n_samples 4 ns = (50 - (26 + 1)) *4 = 148 ns. Is that correct ?

Or should we define some kind of effective window taking into account partial pulses ?

@dneise Can we replace it by :

pulse_mask[:, 1:-1] = (
            (c[:, :-2] <= c[:, 1:-1]) &
            (c[:, 1:-1] >= c[:, 2:]) &
            (c[:, 1:-1] > threshold)
        )

?

If I read this correctly what it does is : verifying

$\\ n \in [0, N_{samples} - 3] \\ x_n<=x_{n+1} \\ x_{n+1} \geq x_{n+2} \\ x_{n+1} \geq T$$

Where the n+1 is where the maximum is right ?

calispac commented 6 years ago

Ah sorry... This has to do with the filter that is applied ?

calispac commented 6 years ago

Maybe one could use this has a smoothing function : https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter1d.html#scipy.ndimage.gaussian_filter1d

dneise commented 6 years ago

from: dark_3MHz_toymc.h5.h5

calispac commented 6 years ago

See this commit https://github.com/cta-sst-1m/digicampipe/pull/182/commits/5279463215014e9ff97766f5e69d4a38132018e4

Where I rewrote the pulse finder with gaussian filter.

dneise commented 6 years ago

Cyril .. I am not sure I can follow you. What is in your opinion the problem?

calispac commented 6 years ago

Cyril .. I am not sure I can follow you. What is in your opinion the problem?

The problem is that right now in find_pulse_3() pulses are found only if they are within [:, 6:-6] therefore we are throwing away some data. But I shouldn't have posted that here since the max-min method is not affected by this. If one evaluates the DCR from the SPE reconstructed with a pulse_finder that means that : DCR = N_pulses / (N_events * T_window) where T_window is therefore shorter

dneise commented 6 years ago

I agree, neither the max_min method, nor the random_trigger method is affected by this.

AndriiNagai commented 6 years ago

Thank you. For max-min, could you do double Gaussian fit in range from 0 to ~12? Because N0 is sitting on the 1 p.e. peak.

Best regards, Andrii

On 11 Apr 2018, at 16:07, Dominik Neise notifications@github.com wrote:

from: dark_3MHz_toymc.h5.h5

https://user-images.githubusercontent.com/8200858/38621963-5caa76d6-3da2-11e8-9c2a-852479ecaae9.png — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-380464826, or mute the thread https://github.com/notifications/unsubscribe-auth/AYRh4BIqgyqSbG3Qp92tr8Rk6pVJWg3_ks5tng4VgaJpZM4TLmcT.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/cta-sst-1m/digicampipe","title":"cta-sst-1m/digicampipe","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/cta-sst-1m/digicampipe"}},"updates":{"snippets":[{"icon":"PERSON","message":"@dneise in #176: from: dark_3MHz_toymc.h5.h5\r\n\r\n\r\n"}],"action":{"name":"View Issue","url":"https://github.com/cta-sst-1m/digicampipe/issues/176#issuecomment-380464826"}}}

cta-sst-1m / digicampipe

Measuring Darkcount rate #176