ECalCalibration/script/timeWalkCor.C

0001 double fit_slewing(double *x, double *par) {
0002   return par[0] + par[1] * pow(x[0] - par[2], par[3]);
0003 }
0004
0005 double fit_landau(double *x, double *par) {
0006   return par[2] * TMath::Landau(x[0], par[0], par[1], kTRUE);
0007 }
0008
0009 void timeWalkCor(TString filename) {
0010
0011   TChain *ch = new TChain("events");
0012   ch->Add(filename);
0013   TTreeReader reader(ch);
0014
0015   // CALOROCHits(p-going)
0016   TTreeReaderArray<int> pos_phase(reader, "EcalBarrelScFiPCALOROCHits.samplePhase");
0017   TTreeReaderArray<int> pos_stamp(reader, "EcalBarrelScFiPCALOROCHits.timeStamp");
0018   TTreeReaderArray<unsigned int> pos_begin(reader, "EcalBarrelScFiPCALOROCHits.bSamples_begin");
0019   TTreeReaderArray<unsigned int> pos_end(reader, "EcalBarrelScFiPCALOROCHits.bSamples_end");
0020   TTreeReaderArray<uint16_t> pos_lowADC( reader, "_EcalBarrelScFiPCALOROCHits_bSamples.lowGainADC");
0021   TTreeReaderArray<uint16_t> pos_TOA(reader, "_EcalBarrelScFiPCALOROCHits_bSamples.timeOfArrival");
0022
0023   // CALOROCHits(e-going)
0024   TTreeReaderArray<int> neg_phase(reader, "EcalBarrelScFiNCALOROCHits.samplePhase");
0025   TTreeReaderArray<int> neg_stamp(reader, "EcalBarrelScFiNCALOROCHits.timeStamp");
0026   TTreeReaderArray<unsigned int> neg_begin(reader, "EcalBarrelScFiNCALOROCHits.bSamples_begin");
0027   TTreeReaderArray<unsigned int> neg_end(reader, "EcalBarrelScFiNCALOROCHits.bSamples_end");
0028   TTreeReaderArray<uint16_t> neg_lowADC(reader, "_EcalBarrelScFiNCALOROCHits_bSamples.lowGainADC");
0029   TTreeReaderArray<uint16_t> neg_TOA(reader, "_EcalBarrelScFiNCALOROCHits_bSamples.timeOfArrival");
0030
0031   // Config variables
0032   const std::size_t n_samples = 7;
0033   const double capADC = 1024;
0034   const double dyRangeADC = 2500;
0035   const double capTOA = 1024;
0036   const double toa_thres = 7;
0037   const double timeWindow = 25;
0038   const double pulse_sigma = 10;
0039
0040   // To perform time-walk correction
0041   TH2D *h2_ADC_dt = new TH2D("h2_ADC_dt", "", 80, 0, 400, 200, -25, 25);
0042   TF1 *f1_slewing = new TF1("f1_slewing", fit_slewing, 5, 500, 4);
0043   f1_slewing->SetNpx(10000);
0044   f1_slewing->SetParameters(-13.7915, 33.5238, 3.15088, -0.313885);
0045
0046   // To reconstruct pulse shapes
0047   TGraph *gr_pulse = new TGraph();
0048   TF1 *f1_landau = new TF1("f1_landau", fit_landau, 0, 300, 3);
0049
0050   // Build the pulse graph, fit it, and fill h2_ADC_dt for one side.
0051   auto processHits = [&](TTreeReaderArray<int> &phase_arr, TTreeReaderArray<int> &stamp_arr,
0052                          TTreeReaderArray<unsigned int> &begin_arr, TTreeReaderArray<unsigned int> &end_arr,
0053                          TTreeReaderArray<uint16_t> &lowADC, TTreeReaderArray<uint16_t> &TOA) {
0054     for (int i = 0; i < phase_arr.GetSize(); i++) {
0055       int begin = begin_arr[i];
0056       int end = end_arr[i];
0057       double phase = phase_arr[i];
0058       double stamp = stamp_arr[i];
0059
0060       int idx_toa = -1;
0061       int idx_toa_rel = -1;
0062       double t_rec = 0;
0063       double adc[n_samples] = {0,};
0064       double adcSum = 0;
0065
0066       for (int j = begin; j < end; j++) {
0067         adcSum += lowADC[j];
0068         adc[j - begin] = lowADC[j];
0069
0070         if (TOA[j] > 0) {
0071           idx_toa = j;
0072           idx_toa_rel = j - begin;
0073         }
0074       }
0075
0076       if (idx_toa > -1 && adcSum > toa_thres) {
0077         // Add (t_rec, toa_thres) point from TOA
0078         t_rec = (phase - TOA[idx_toa]) * (timeWindow / capTOA) +
0079                 (stamp + idx_toa_rel) * timeWindow;
0080         gr_pulse->Set(0);
0081         gr_pulse->SetPoint(0, t_rec, toa_thres);
0082
0083         // Add (t, Npe) points from ADC samples
0084         for (int j = idx_toa_rel; j < n_samples; j++) {
0085           if (adc[j] > 0)
0086             gr_pulse->SetPoint(gr_pulse->GetN(),
0087                                t_rec + TOA[idx_toa] * (timeWindow / capTOA) +
0088                                    timeWindow * (j - idx_toa_rel),
0089                                adc[j] * (dyRangeADC / capADC));
0090         }
0091
0092         // Fit pulse shape
0093         f1_landau->SetParameters(t_rec + 20, pulse_sigma, adcSum * 55.37);
0094         f1_landau->SetParLimits(1, pulse_sigma - 2, pulse_sigma + 2);
0095         int status = gr_pulse->Fit("f1_landau", "RBQ0");
0096         if (status < 0)
0097           continue;
0098
0099         double chi2 = f1_landau->GetChisquare();
0100         double npe_fit = f1_landau->GetMaximum();
0101         double t_fit = f1_landau->GetX(0.3 * npe_fit, 0, f1_landau->GetParameter(0));
0102         if (std::isnan(t_fit))
0103           t_fit = 0;
0104
0105         if (adcSum > 0 && t_rec > 0 && t_fit > 0 && chi2 < 2) {
0106           h2_ADC_dt->Fill(adcSum, t_rec - t_fit);
0107         }
0108       }
0109     }
0110   };
0111
0112   int nevent = 0;
0113   while (reader.Next()) {
0114     nevent++;
0115     if (nevent % 500 == 0)
0116       printf(">>> %d\n", nevent);
0117
0118     processHits(pos_phase, pos_stamp, pos_begin, pos_end, pos_lowADC, pos_TOA);
0119     processHits(neg_phase, neg_stamp, neg_begin, neg_end, neg_lowADC, neg_TOA);
0120   }
0121
0122   TCanvas *canv = new TCanvas("canv", "", 600, 600);
0123   h2_ADC_dt->Draw();
0124   h2_ADC_dt->Fit("f1_slewing", "R");
0125 }