ECalCalibration/script/zRec.C

0001 void zRec(TString filename) {
0002
0003   TChain *ch = new TChain("events");
0004   ch->Add(filename);
0005   TTreeReader reader(ch);
0006
0007   // True z position
0008   TTreeReaderArray<uint64_t> id(reader, "EcalBarrelScFiPAttenuatedHits.cellID");
0009   TTreeReaderArray<float> z_true(reader, "EcalBarrelScFiPAttenuatedHits.position.z");
0010
0011   // CALOROCHits(p-going)
0012   TTreeReaderArray<uint64_t> pos_id(reader, "EcalBarrelScFiPCALOROCHits.cellID");
0013   TTreeReaderArray<int> pos_phase(reader, "EcalBarrelScFiPCALOROCHits.samplePhase");
0014   TTreeReaderArray<int> pos_stamp(reader, "EcalBarrelScFiPCALOROCHits.timeStamp");
0015   TTreeReaderArray<unsigned int> pos_begin(reader, "EcalBarrelScFiPCALOROCHits.bSamples_begin");
0016   TTreeReaderArray<unsigned int> pos_end(reader, "EcalBarrelScFiPCALOROCHits.bSamples_end");
0017   TTreeReaderArray<uint16_t> pos_lowADC(reader, "_EcalBarrelScFiPCALOROCHits_bSamples.lowGainADC");
0018   TTreeReaderArray<uint16_t> pos_TOA(reader, "_EcalBarrelScFiPCALOROCHits_bSamples.timeOfArrival");
0019
0020   // CALOROCHits(e-going)
0021   TTreeReaderArray<uint64_t> neg_id(reader, "EcalBarrelScFiNCALOROCHits.cellID");
0022   TTreeReaderArray<int> neg_phase(reader, "EcalBarrelScFiNCALOROCHits.samplePhase");
0023   TTreeReaderArray<int> neg_stamp(reader, "EcalBarrelScFiNCALOROCHits.timeStamp");
0024   TTreeReaderArray<unsigned int> neg_begin(reader, "EcalBarrelScFiNCALOROCHits.bSamples_begin");
0025   TTreeReaderArray<unsigned int> neg_end(reader, "EcalBarrelScFiNCALOROCHits.bSamples_end");
0026   TTreeReaderArray<uint16_t> neg_lowADC(reader, "_EcalBarrelScFiNCALOROCHits_bSamples.lowGainADC");
0027   TTreeReaderArray<uint16_t> neg_TOA(reader, "_EcalBarrelScFiNCALOROCHits_bSamples.timeOfArrival");
0028
0029   // Config variables
0030   const double capTOA = 1024;
0031   const double toa_thres = 7;
0032   const double timeWindow = 25;
0033
0034   // z position of the e-going end
0035   const double z0_neg = -2637.5;
0036   // Time-walk correction parameters from timeWalkCor.C
0037   const double pars[4] = {-13.7915, 33.5238, 3.15088, -0.313885};
0038
0039   // To extract z_true
0040   TH2D* h2_dt_z = new TH2D("h2_dt_z", "", 200, -20, 30, 100, 500, 4500);
0041   TF1* f1_dt_z = new TF1("f1_dt_z", "[0]*x+[1]", -15, 25);
0042   f1_dt_z->SetParameters(83.3221, 2219.58);
0043   std::unordered_map<uint64_t, double> id_t_pos, id_t_neg, id_z_true;
0044
0045   // Reconstruct the time-walk-corrected time of each hit and store it by cellID.
0046   auto processHits = [&](std::unordered_map<uint64_t, double> &id_t,
0047                  TTreeReaderArray<uint64_t> &id_arr,
0048              TTreeReaderArray<int> &phase_arr, TTreeReaderArray<int> &stamp_arr,
0049                  TTreeReaderArray<unsigned int> &begin_arr, TTreeReaderArray<unsigned int> &end_arr,
0050                  TTreeReaderArray<uint16_t> &lowADC, TTreeReaderArray<uint16_t> &TOA) {
0051     for (int i = 0; i < id_arr.GetSize(); i++) {
0052       uint64_t cellID = id_arr[i];
0053       int begin = begin_arr[i];
0054       int end = end_arr[i];
0055       double phase = phase_arr[i];
0056       double stamp = stamp_arr[i];
0057
0058       int idx_toa = -1;
0059       int idx_toa_rel = -1;
0060       double adcSum = 0;
0061
0062       for (int j = begin; j < end; j++) {
0063         adcSum += lowADC[j];
0064
0065         if (TOA[j] > 0) {
0066           idx_toa = j;
0067           idx_toa_rel = j - begin;
0068         }
0069       }
0070
0071       if (idx_toa > -1 && adcSum > toa_thres) {
0072         double t_rec = (phase - TOA[idx_toa]) * (timeWindow / capTOA) +
0073                        (stamp + idx_toa_rel) * timeWindow;
0074         id_t[cellID] = t_rec - (pars[1] * pow(adcSum - pars[2], pars[3]) + pars[0]);
0075       }
0076     }
0077   };
0078
0079   int nevent = 0;
0080   while (reader.Next()) {
0081     nevent++;
0082     if (nevent % 500 == 0)
0083       printf(">>> %d\n", nevent);
0084
0085     id_t_pos.clear();
0086     id_t_neg.clear();
0087     id_z_true.clear();
0088
0089     for (int i = 0; i < id.GetSize(); i++)
0090       id_z_true[id[i]] = z_true[i];
0091
0092     processHits(id_t_pos, pos_id, pos_phase, pos_stamp, pos_begin, pos_end,
0093                 pos_lowADC, pos_TOA);
0094     processHits(id_t_neg, neg_id, neg_phase, neg_stamp, neg_begin, neg_end,
0095                 neg_lowADC, neg_TOA);
0096
0097     for (auto &[cellID, t_pos] : id_t_pos) {
0098       auto it = id_t_neg.find(cellID);
0099       if (it != id_t_neg.end())
0100         h2_dt_z->Fill(it->second - t_pos, id_z_true[cellID] - z0_neg);
0101     }
0102   }
0103
0104   TCanvas *canv = new TCanvas("canv", "", 600, 600);
0105   h2_dt_z->Draw();
0106   h2_dt_z->Fit("f1_dt_z", "R");
0107 }