ECalCalibration/script/noisePars.C

0001 #include <DDDigi/noise/FalphaNoise.h>
0002
0003 double pulse(double *x, double *par) {
0004   return par[0] * TMath::Landau(x[0], par[2] * par[1], par[1], kTRUE);
0005 }
0006
0007 void setGraphStyle(TGraph* gr, Style_t style,  Color_t color, Size_t size);
0008
0009 void noisePars() {
0010
0011   // Signal
0012   const double npe_mean = 100;
0013   const double npe_sigma = 5;
0014
0015   // To draw pulses for comparison
0016   constexpr double t_min = -300;
0017   constexpr double t_max = 300;
0018   constexpr double dt = 0.5;
0019   constexpr int n_bins = static_cast<int>((t_max - t_min) / dt);
0020
0021   // Config variables
0022   const double gain = 55.37;
0023   const double sigma_offset = 3.5;
0024   const double sigma_analog = 10;
0025
0026   const int n_samples = 7;
0027   const double phase = 20;
0028   const double time_window = 25;
0029
0030   // Histograms
0031   TH1D* h1_adc_real = new TH1D("h1_adc_real", "", 100, 0, 450);
0032   TH1D* h1_adc_model = new TH1D("h1_adc_model", "", 100, 0, 450);
0033
0034   // Real noise
0035   TRandom3 rng(0);
0036   const double dark_rate = 112.97e6;
0037   const double window_sec = (t_max - t_min) * 1.0e-9;
0038   const double mean_n_dark = dark_rate * window_sec;
0039
0040   // Falpha noise
0041   const int poles = 5;
0042   const double variance = 1.0;
0043   const double alpha = 1.8;
0044   const double scale = 0.35;
0045   const double offset = 6.1;
0046   std::default_random_engine gen(12345);
0047   dd4hep::detail::FalphaNoise falpha(poles, alpha, variance);
0048
0049   const int nev = 10000;
0050
0051   // Example pulse shapes at a specific event
0052   const int evnum = 30;
0053   double times[n_bins] = {0};
0054   double amps_real[n_bins] = {0};
0055   double amps_model[n_bins] = {0};
0056   double sig_amps[n_bins] = {0};
0057   double noise_amps_real[n_bins] = {0};
0058   double noise_amps_model[n_bins] = {0};
0059
0060   for (int i = 0; i < nev; i++) {
0061     if ((i % 500) == 0)
0062       cout << i << endl;
0063
0064     const double sig_npe = std::max(0.0, rng.Gaus(npe_mean, npe_sigma));
0065     double adc_real = 0;
0066     double adc_model = 0;
0067
0068     const int n_dark = rng.Poisson(mean_n_dark);
0069     std::vector<double> dark_times;
0070     std::vector<double> dark_amps;
0071
0072     for (int i = 0; i < n_dark; i++) {
0073       // Dark noises are generated n_dark times randomly.
0074       dark_times.push_back(rng.Uniform(t_min, t_max));
0075       // It has 1 p.e. pulse height.
0076       dark_amps.push_back(1);
0077     }
0078
0079     for (int j = 0; j < n_bins; j++) {
0080       double time = t_min + j * dt;
0081       double sig = 0;
0082       sig += (sig_npe * gain) * TMath::Landau(time, sigma_analog * sigma_offset,
0083                                               sigma_analog, kTRUE);
0084       double noise_real = 0;
0085       for (int k = 0; k < n_dark; k++) {
0086         noise_real +=
0087             (dark_amps[k] * gain) *
0088             TMath::Landau(time, dark_times[k] + sigma_analog * sigma_offset,
0089                           sigma_analog, kTRUE);
0090       }
0091
0092       double noise_model = std::max(0.0, scale * falpha(gen) + offset);
0093
0094       // Store pulse shapes for comparison
0095       if (i == evnum) {
0096         times[j] = time;
0097         sig_amps[j] = sig;
0098         noise_amps_real[j] = noise_real;
0099         noise_amps_model[j] = noise_model;
0100         amps_real[j] = sig + noise_real;
0101         amps_model[j] = sig + noise_model;
0102       }
0103
0104       // CALOROC measurement
0105       for (int k = 0; k < n_samples; k++) {
0106         if (time == phase + k * time_window) {
0107             adc_real += sig + noise_real;
0108             adc_model += sig + noise_model;
0109         }
0110       }
0111     }
0112     h1_adc_real->Fill(adc_real);
0113     h1_adc_model->Fill(adc_model);
0114   }
0115
0116   // For example pulses comparison
0117   TH2D *h2_pulse = new TH2D("h2_pulse", "", 6, -100, t_max, 6, 0, npe_mean + 4 * npe_sigma);
0118
0119   TGraph *gr_sig = new TGraph(n_bins, times, sig_amps);
0120   setGraphStyle(gr_sig, 20, 4, 0.3);
0121
0122   TGraph *gr_noise_real = new TGraph(n_bins, times, noise_amps_real);
0123   setGraphStyle(gr_noise_real, 20, 2, 0.3);
0124   TGraph *gr_noise_model = new TGraph(n_bins, times, noise_amps_model);
0125   setGraphStyle(gr_noise_model, 20, 2, 0.3);
0126
0127   TGraph *gr_pulse_real = new TGraph(n_bins, times, amps_real);
0128   setGraphStyle(gr_pulse_real, 20, 1, 0.3);
0129   TGraph *gr_pulse_model = new TGraph(n_bins, times, amps_model);
0130   setGraphStyle(gr_pulse_model, 20, 1, 0.3);
0131
0132   // Example pulse shapes
0133   TCanvas *c_pulse = new TCanvas("c_pulse", "", 1300, 1000);
0134   c_pulse->Divide(1, 2);
0135
0136   c_pulse->cd(1);
0137   h2_pulse->Draw();
0138   gr_sig->Draw("PL");
0139   gr_noise_real->Draw("PL");
0140   gr_pulse_real->Draw("PL");
0141
0142   c_pulse->cd(2);
0143   h2_pulse->Draw();
0144   gr_sig->Draw("PL");
0145   gr_noise_model->Draw("PL");
0146   gr_pulse_model->Draw("PL");
0147
0148   // ADC distribution comparison
0149   TCanvas *c_adc = new TCanvas("c_adc", "", 1200, 600);
0150   c_adc->Divide(2, 1);
0151
0152   c_adc->cd(1);
0153   h1_adc_real->Draw();
0154
0155   c_adc->cd(2);
0156   h1_adc_model->Draw();
0157 }
0158
0159 void setGraphStyle(TGraph* gr, Style_t style,  Color_t color, Size_t size){
0160
0161   gr->SetMarkerStyle(style);
0162   gr->SetMarkerSize(size);
0163   gr->SetMarkerColor(color);
0164   gr->SetLineColor(color);
0165 }