mbd/TIMING/make_mbd_l1lut.C

0001 //
0002 // Make the MBD LUT's
0003 //
0004 #include <mbd/MbdCalib.h>
0005 #include <mbd/MbdDefs.h>
0006 #include <mbd/MbdGeomV1.h>
0007
0008 #include <Rtypes.h>  // defines R__LOAD_LIBRARY macro for clang-tidy
0009 #include <TGraph.h>
0010 #include <TMath.h>
0011
0012 #include <cmath>
0013 #include <cstdint>
0014 #include <fstream>
0015
0016 R__LOAD_LIBRARY(libmbd.so)
0017
0018 MbdGeom *mbdgeom{nullptr};
0019 MbdCalib *mcal{nullptr};
0020
0021 float tdclut[128][1024] = {};
0022 uint16_t adclut[128][1024] = {};
0023 uint16_t slewlut[128][4096] = {};
0024 double step = 0;  // ns/LL1 tdc count
0025 double slewcorr_step[128] = {};
0026
0027 void make_generic_adclut();
0028 void make_generic_slewlut();
0029 void make_adclut(const char *scorr_fname = "mbd_slewcorr.calib");
0030 void make_slewlut(const char *scorr_fname = "mbd_slewcorr.calib");
0031
0032 void make_generic_adclut()
0033 {
0034   // write adc lut
0035   std::ofstream adclutfile("mbdadc.lut");
0036   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0037   {
0038     adclutfile << "adc " << ipmt << std::endl;
0039     for (unsigned int iaddr = 0; iaddr < 1024; iaddr++)
0040     {
0041       uint16_t qsum7 = iaddr >> 3U;
0042       uint16_t slew3 = iaddr >> 7U;
0043       adclut[ipmt][iaddr] = (slew3 << 7U) + qsum7;  // assuming slew3 are high order bits
0044
0045       adclutfile << slew3 * 1000 + qsum7 << "\t";
0046
0047       if (iaddr % 8 == 7)
0048       {
0049         adclutfile << std::endl;
0050       }
0051     }
0052   }
0053   adclutfile.close();
0054 }
0055
0056 void make_generic_slewlut()
0057 {
0058   // write slew lut
0059   std::ofstream slewlutfile("mbdslew.lut");
0060   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0061   {
0062     slewlutfile << "slew " << ipmt << std::endl;
0063     for (int islew = 0; islew < 8; islew++)  // 3-bit
0064     {
0065       //      uint16_t saddr = islew;
0066       for (int itdc = 0; itdc < 512; itdc++)  // 9-bit
0067       {
0068         uint16_t taddr = itdc;
0069         if (taddr != 0)
0070         {
0071           slewlut[ipmt][itdc] = 1000 + taddr;
0072           // 1000 represents the active bit in the lut file
0073           // when loading into the actual MBD ADC LUT, we should
0074           // we should set bit 10 = 1.
0075         }
0076         else
0077         {
0078           slewlut[ipmt][itdc] = taddr;
0079         }
0080
0081         slewlutfile << slewlut[ipmt][itdc] << "\t";
0082
0083         if (itdc % 8 == 7)
0084         {
0085           slewlutfile << std::endl;
0086         }
0087       }
0088     }
0089   }
0090   slewlutfile.close();
0091 }
0092
0093 void make_adclut(const char * /*scorr_fname*/)
0094 {
0095   // write adc lut
0096   std::ofstream adclutfile("mbdadc.lut");
0097   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0098   {
0099     int feech = mbdgeom->get_feech(ipmt, 0);
0100
0101     // get the max slew correction, which we set at ADC=40
0102     double max_slew = mcal->get_scorr(feech, 40);
0103     slewcorr_step[ipmt] = max_slew / 8;
0104     //    int nsteps = 0;
0105
0106     adclutfile << "adc " << ipmt << std::endl;
0107     for (unsigned int iaddr = 0; iaddr < 1024; iaddr++)
0108     {
0109       uint16_t qsum7 = iaddr >> 3U;
0110
0111       uint16_t actual_adc = (iaddr << 4U) + 8;  // 8 = 2^4/2, so actual_adc is at midpoint
0112       double scorr = mcal->get_scorr(feech, actual_adc);
0113       int scorr_step = static_cast<int>(scorr / slewcorr_step[ipmt]);
0114       if (scorr_step < 0)
0115       {
0116         scorr_step = 0;
0117       }
0118       else if (scorr_step > 7)
0119       {
0120         scorr_step = 7;
0121       }
0122
0123       uint16_t slew3 = 7 - scorr_step;
0124
0125       if (ipmt == 0)
0126       {
0127         std::cout << "adcslewbin " << iaddr << "\t" << scorr << "\t" << scorr_step << "\t" << slew3 << std::endl;
0128       }
0129
0130       adclut[ipmt][iaddr] = (slew3 << 7U) + qsum7;  // assuming slew3 are high order bits
0131
0132       adclutfile << slew3 * 1000 + qsum7 << "\t";
0133
0134       if (iaddr % 8 == 7)
0135       {
0136         adclutfile << std::endl;
0137       }
0138     }
0139   }
0140   adclutfile.close();
0141 }
0142
0143 void make_slewlut(const char * /*scorr_fname*/)
0144 {
0145   // write slew lut
0146   std::ofstream slewlutfile("mbdslew.lut");
0147   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0148   {
0149     slewlutfile << "slew " << ipmt << std::endl;
0150     for (int islew = 0; islew < 8; islew++)  // 3-bit
0151     {
0152       //      uint16_t saddr = islew;
0153       for (int itdc = 0; itdc < 512; itdc++)  // 9-bit
0154       {
0155         int taddr = itdc;
0156         if (taddr != 0)
0157         {
0158           // apply slew correction
0159           taddr -= static_cast<int>((7 - islew) * (slewcorr_step[ipmt] / step));
0160           if (ipmt == 0)
0161           {
0162             std::cout << "slewxxx " << islew << "\t" << itdc << "\t" << (7 - islew) << "\t" << slewcorr_step[ipmt] << "\t" << step << "\t" << (slewcorr_step[ipmt] / step) << "\t" << taddr << std::endl;
0163           }
0164
0165           if (taddr < 0)
0166           {
0167             taddr = 1;
0168           }
0169           slewlut[ipmt][itdc] = 1000 + taddr;
0170           // 1000 represents the active bit in the lut file
0171           // when loading into the actual MBD ADC LUT, we should
0172           // we should set bit 10 = 1.
0173         }
0174         else  // 0 is no-hit for TDC
0175         {
0176           slewlut[ipmt][itdc] = taddr;
0177         }
0178
0179         slewlutfile << slewlut[ipmt][itdc] << "\t";
0180
0181         if (itdc % 8 == 7)
0182         {
0183           slewlutfile << std::endl;
0184         }
0185       }
0186     }
0187   }
0188   slewlutfile.close();
0189 }
0190
0191 void make_mbd_l1lut(const char *tcorr_fname = "mbd_timecorr.calib",
0192                     const char *t0_fname = "mbd_tt_t0.calib",
0193                     const char *scorr_fname = "mbd_slewcorr.calib",
0194                     const char *tpranges_fname = "mbd_tpscan.txt")
0195 {
0196   mbdgeom = new MbdGeomV1();
0197
0198   mcal = new MbdCalib();
0199   mcal->Download_TimeCorr(tcorr_fname);
0200   mcal->Download_TTT0(t0_fname);
0201   mcal->Download_SlewCorr(scorr_fname);
0202
0203   std::ifstream tprangesfile(tpranges_fname);
0204   int temp_feech;
0205   double mintdc[128] = {0};
0206   double maxtdc[128] = {0};
0207   double mintime[128] = {0};
0208   double maxtime[128] = {0};
0209   double dtdc[128] = {0};
0210   double dtime[128] = {0};
0211
0212   // Find the overlapping good time ranges for all the channels
0213   TGraph *g_mintdc = new TGraph();
0214   TGraph *g_maxtdc = new TGraph();
0215   g_mintdc->SetName("g_mintdc");
0216   g_maxtdc->SetName("g_maxtdc");
0217   TGraph *g_mintimes = new TGraph();
0218   TGraph *g_maxtimes = new TGraph();
0219   g_mintimes->SetName("g_mintdc");
0220   g_maxtimes->SetName("g_maxtdc");
0221
0222   while (tprangesfile >> temp_feech)
0223   {
0224     int pmtch = (double) mbdgeom->get_pmt(temp_feech);
0225
0226     tprangesfile >> mintdc[pmtch] >> maxtdc[pmtch] >> dtdc[pmtch] >> mintime[pmtch] >> maxtime[pmtch] >> dtime[pmtch];
0227
0228     // std::cout << "maxtime " << pmtch << "\t" << maxtime[pmtch] << "\t" << mcal->get_tcorr( temp_feech, 16 ) << std::endl;;
0229     maxtime[pmtch] = mcal->get_tcorr(temp_feech, 16) - mcal->get_tt0(pmtch);
0230
0231     std::cout << "mintime " << pmtch << "\t" << mintime[pmtch] << "\t" << mcal->get_tcorr(temp_feech, static_cast<int>(std::round(maxtdc[pmtch]))) << std::endl;
0232     ;
0233     std::cout << mcal->get_tt0(pmtch) << std::endl;
0234     mintime[pmtch] = mcal->get_tcorr(temp_feech, static_cast<int>(round(maxtdc[pmtch]))) - mcal->get_tt0(pmtch);
0235
0236     // correct for bad channels
0237     if (std::abs(mcal->get_tt0(pmtch)) > 100.)
0238     {
0239       maxtime[pmtch] = mcal->get_tcorr(temp_feech, 16) - 3.0;
0240       mintime[pmtch] = mcal->get_tcorr(temp_feech, static_cast<int>(round(maxtdc[pmtch]))) - 3.0;
0241     }
0242
0243     int n = g_mintdc->GetN();
0244     g_mintdc->SetPoint(n, pmtch, mintdc[pmtch]);
0245
0246     n = g_maxtdc->GetN();
0247     g_maxtdc->SetPoint(n, pmtch, maxtdc[pmtch]);
0248
0249     n = g_mintimes->GetN();
0250     g_mintimes->SetPoint(n, pmtch, mintime[pmtch]);
0251
0252     n = g_maxtimes->GetN();
0253     g_maxtimes->SetPoint(n, pmtch, maxtime[pmtch]);
0254   }
0255
0256   int n = g_mintimes->GetN();
0257   double lut_mintime = TMath::MaxElement(n, g_mintimes->GetY());
0258
0259   n = g_maxtimes->GetN();
0260   double lut_maxtime = TMath::MinElement(n, g_maxtimes->GetY());
0261
0262   std::cout << "max-min time: " << lut_mintime << "\t" << lut_maxtime << "\t" << lut_maxtime - lut_mintime << std::endl;
0263
0264   g_mintimes->SetMarkerStyle(20);
0265   g_mintimes->SetMarkerColor(3);
0266   g_maxtimes->SetMarkerStyle(20);
0267   g_maxtimes->SetMarkerColor(4);
0268   g_mintimes->SetMinimum(lut_mintime * 1.5);
0269   g_mintimes->SetMaximum(lut_maxtime * 1.2);
0270   g_mintimes->Draw("ap");
0271   g_maxtimes->Draw("p");
0272
0273   double max_range = lut_maxtime - lut_mintime;
0274   step = max_range / 511.;  // ns/LL1 tdc count
0275   double zstep = (2.0 / 29.9792458) * (1.0 / step);
0276   std::cout << "L1 trig time unit is " << step << " ns/count" << std::endl;
0277   std::cout << "L1 trig time unit is " << zstep << " count/cm" << std::endl;
0278
0279   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0280   {
0281     if (mbdgeom->get_type(ifeech) == 1)
0282     {
0283       continue;  // skip q-ch's
0284     }
0285     int pmtch = mbdgeom->get_pmt(ifeech);
0286
0287     for (unsigned int iaddr = 0; iaddr < 1024; iaddr++)
0288     {
0289       // could average over the 16 values
0290       unsigned int tdc = (iaddr << 4U) + 8;
0291       float true_time = mcal->get_tcorr(ifeech, tdc) - mcal->get_tt0(pmtch) - lut_mintime;
0292
0293       if (mcal->get_tt0(pmtch) < -100.)
0294       {
0295         // set bad channels to 0;
0296         // tdclut[pmtch][iaddr] = 0;
0297         // continue;
0298
0299         true_time = mcal->get_tcorr(ifeech, tdc) - 3.0 - lut_mintime;
0300       }
0301
0302       if (true_time < 0. || true_time > max_range || std::isnan(true_time) || true_time > 22.5)
0303       {
0304         tdclut[pmtch][iaddr] = 0;
0305       }
0306       else
0307       {
0308         if (std::isnan(true_time))
0309         {
0310           std::cout << "std::isnan " << ifeech << " " << mcal->get_tcorr(ifeech, tdc)
0311                     << "\t" << mcal->get_tt0(pmtch) << std::endl;
0312           ;
0313         }
0314
0315         // tdclut[pmtch][iaddr] = static_cast<int>( round(true_time/step) );
0316         tdclut[pmtch][iaddr] = round(true_time / step);
0317       }
0318
0319       if (ifeech == 0 && (iaddr < 24 || iaddr > (1023 - 96)))
0320       {
0321         if (iaddr % 8 == 0)
0322         {
0323           std::cout << iaddr << "\t";
0324         }
0325         std::cout << tdclut[pmtch][iaddr] << "\t";  // be sure to set upper bit to 1 for non-zero values when writing to actual LUT in ADC
0326         if (iaddr % 8 == 7)
0327         {
0328           std::cout << std::endl;
0329         }
0330       }
0331     }
0332   }
0333
0334   // write tdc lut
0335   std::ofstream tdclutfile("mbdtdc.lut");
0336   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0337   {
0338     if (mbdgeom->get_type(ifeech) == 1)
0339     {
0340       continue;  // skip q-ch's
0341     }
0342     int pmtch = mbdgeom->get_pmt(ifeech);
0343
0344     tdclutfile << "tdc " << pmtch << std::endl;
0345     for (int iaddr = 0; iaddr < 1024; iaddr++)
0346     {
0347       tdclutfile << tdclut[pmtch][iaddr] << "\t";
0348
0349       if (iaddr % 8 == 7)
0350       {
0351         tdclutfile << std::endl;
0352       }
0353     }
0354   }
0355   tdclutfile.close();
0356
0357   make_adclut();
0358   make_slewlut();
0359 }